UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Solving reachable sets on a manifold Cross, Elizabeth Ann 2007

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.

Item Metadata

Download

Media
831-ubc_2007-0371.pdf [ 6.48MB ]
Metadata
JSON: 831-1.0052081.json
JSON-LD: 831-1.0052081-ld.json
RDF/XML (Pretty): 831-1.0052081-rdf.xml
RDF/JSON: 831-1.0052081-rdf.json
Turtle: 831-1.0052081-turtle.txt
N-Triples: 831-1.0052081-rdf-ntriples.txt
Original Record: 831-1.0052081-source.json
Full Text
831-1.0052081-fulltext.txt
Citation
831-1.0052081.ris

Full Text

r^-Solving Reachable Sets on a Manifold by Elizabeth Ann Cross B.Math., Carleton University, 2005 A THESIS S U B M I T T E D IN P A R T I A L F U L F I L M E N T OF T H E R E Q U I R E M E N T S F O R T H E D E G R E E OF Master of Science in The Faculty of Graduate Studies (Computer Science) The University Of British Columbia August 30, 2007 © Elizabeth Ann Cross 2007 11 Abstract First, this thesis explores the implementation of the fast marching method as part of the toolbox of level set methods. This method uses Dijkstra's algo-rithm to approximate the solution to the non-linear Eikonal equation. Functions for calculating signed distances and extension velocities are also implemented. These functions use the fast marching method in their implementation. Second, it explores a method for computing reachable sets on a manifold; in other words, the dynamics governing these reachable sets can be described by a Differential Algebraic Equation. It uses level set methods to solve the underlying Hamilton Jacobi equation of the reachable set and it ensures an accurate solution on the manifold by using the closest point method. The closest point method guarantees that the reachable set is perpendicular to the manifold at the points of intersection. Several two and three dimensional toy problems and a real-life power generator problem are explored in order to test the method. i i i Contents Abstract i i Contents i i i List of Tables v List of Figures vi Acknowledgements vii 1 Introduction 1 1.1 Implicit Surface Functions . . . . ' 1 2 Fast Marching Methods 3 2.1 Background 3 2.1.1 Fast Marching Method 3 2.1.2 Signed Distance Function 5 2.1.3 Extension Velocity 8 2.2 Implementation 9 2.2.1 Implementation Choices 10 2.2.2 Dijkstra • 10 2.2.3 Signed Distance 10 2.2.4 Extension Velocity and Signed Distance 11 2.3 Examples 11 2.3.1 Shape From Shading 11 2.3.2 Extension Velocity of a Circle 13 2.3.3 Movement of a Circle in Normal Direction 15 3 Reachable Sets on a Manifold 18 3.1 Background 18 3.1.1 Differential Algebraic Equations 18 3.1.2 Closest Point Method 18 3.1.3 Reachable Sets 19 3.2 Implementation 21 3.2.1 Level Set Methods 21 3.2.2 Differential Algebraic Equations 21 3.2.3 Closest Point 21 \ Contents iv 3.3 Examples 24 3.3.1 Two-dimensional D A E Toy Problem 24 3.3.2 Three-dimensional D A E Toy Problem with Codimension 1 25 4 Predicting Voltage Instability of a Power System 34 4.1 Background 34 4.2 Closest Point Method used as an Extension Velocity 35 4.3 Implementation 38 4.4 Examples 38 5 Conclusion 48 Bibliography . 5 0 V List of Tables 3.1 Description of parameters used in the closest point operator al-gorithm 22 3.2 Description of functions and their input parameters used in the closest point operator algorithm 23 3.3 Description of functions and their output parameters used in the closest point operator, algorithm 23 3.4 Interpolated implicit surface value for each timestep of the two-dimensional toy D A E problem 26 3.5 Convergence of error for the two dimensional toy example . . . . 26 3.6 The average interpolated implicit surface value for multiple O D E trajectories at each timestep of the three-dimensional toy D A E problem 28 3.7 The maximum interpolated implicit surface value for multiple O D E trajectories at each timestep of the three-dimensional toy D A E problem : 28 4.1 Physical Meaning of Variables and Parameters in Nonlinear Hy-brid Automaton 36 V List of Figures 1.1 Implicit surface representation of an arbitrary curve 2 2.1 Dijkstra's method on a undirected discrete graph 4 2.2 Intersection of the interface with grid lines 7 2.3 Smooth solution of the shape from shading problem 12 2.4 Continuous but not smooth solution of the shape from shading problem 13 2.5 Original and extension velocity functions of a circle . • 14 2.6 Initial and signed distance functions of a circle 14 2.7 Contour plots for a circle moving in its normal direction 16 2.8 Ordered gradient magnitude plots 17 3.1 Reachable set plots for the two dimensional toy example 29 3.2 Convergence rates for the two dimensional toy example 30 3.3 Reachable set plots for the three dimensional toy example in 3D 31 3.4 Plot of the reachable set on the manifold in x2 versus £3 dimen-sions for the three dimensional toy example 32 3.5 Plot of the reachable set on the manifold in s, the horizontal distance along the manifold from the origin in the x2 dimension, versus X3 dimensions for three dimensional toy example 33 4.1 Target set and singular set plots for the power generator problem 40 4.2 . Trajectory simulations generated in three different ways by Mat-lab's D A E / O D E solver 41 4.3 Sample trajectories using the interpolated extension velocity . . . 42 4.4 View of reach tube at t = 0 and t = 1 43 4.5 View of reach tube at t = 2 and t = 3 44 4.6 View of reach tube at t = 4 and t = 5 45 4.7 Several views of the reach tube at the final timestep 46 4.8 The plot compares the interface for the reachable tube at the fi-nal timestep with simulations calculated by Matlab's D A E solver using the dynamics described in (4.1) 47 V l l Acknowledgements First and foremost I would like to thank Ian Mitchell for providing me with excellent guidance and support throughout my thesis. He also kept me organized and on track during the entire thesis process. I would also like to thank Robert Bridson for providing insightful comments about my thesis. Next I would like to thank my SCL labmates for scientific and not so scientific discussions. These discussions made everyday in the lab an enjoyable day. I would also like to thank my other friends in the department and my friends at U B C . The courses I took in my first year provided me with a wide variety of scientific computing knowledge. I would like to thank the professors who taught these interesting courses. Finally I would like to thank my moral support system: Heather, Doug, Angie and Clint. They have always supported me through whatever new ad-ventures I chose to pursue. 1 Chapter 1 Introduction The first purpose of this research is to add some functionality to Ian Mitchell's Level Set Toolbox [6], and this functionality is an implementation of the fast marching method. This method uses Dijkstra's algorithm to approximate the solution to the non-linear Eikonal equation. In addition, functions for calcu-lating signed distances and extension velocities are also implemented. The fast marching method is used in the calculation of both signed distances and exten-sion velocities. The second purpose of this research is to develop a method for solving reach-able sets on a manifold. The dynamics governing the reachable set evolution can be described by a Differential Algebraic Equation (DAE), where the manifold is an algebraic constraint and the movement of the reachable set is described by a differential constraint. A common method used to solve reachable sets in a full dimensional space is the level set method. In this case the reachable set is represented by an implicit surface function. This research adapts these level set techniques in order to solve reachable sets on a manifold. This adapta-tion uses the closest point method to ensure that the implicit surface function representing the reachable set is perpendicular to the manifold at all points of intersection. The function for calculating signed distances is used in the closest point method implementation. Thus, the work in the first part of the thesis is incorporated into second part. This thesis will describe fast marching methods, signed distance functions and extension velocities. It will then outline the implementation of these meth-ods and show some example problems to demonstrate that the implementation is correct. Next, this thesis will outline the method developed to solve reach-able sets on a manifold. It will also demonstrate the method by solving two toy problems. Finally, the thesis will describe a power generator example and demonstrate how this real life D A E can be solved with the previously described methods. 1.1 Implicit Surface Functions In this thesis, closed surfaces are represented by implicit surface functions, <j> [7]. A n interface with a dimension of n — 1 is embedded as an isocontour into a function cj> with an input dimension of n and an output dimension of 1. The zero isocontour, <f>{x) = 0, is usually used to represent the interface but any isocontour, <p{x) = a, can be used for this purpose. These isocontours are equal Chapter 1. Introduction 2 r n + (j) = 0 0>O interface outside Figure 1.1: Implicit surface representation of an arbitrary curve. up to a scalar translation, <f>(x) = cj>(x) — a. This thesis uses the zero isocontour to represent interfaces, so the explicit representation of the interface is defined as T = {x | 4>(x) = 0}. The zero isocontour divides the domain into two sets. The set inside the interface is defined by Q,~ = {x \ cj)(x) < 0} and the set outside the interface is defined by f2 + = {x | <f>(x) > 0}. The majority of the typical geometric tasks for sets and surfaces are easily done with the implicit surface representation. This representation makes it simple to identify which side of the interface a point is on simply by looking at the sign of 4>(x). It is also easy to construct representations by combining several implicit surface representations together. For example, the intersection of the interiors of two interfaces can be represented by taking the maximum value between the two implicit surface representations at each point in the domain, 4>z{x) = m&x.(4>\{x), <f>2{x)). The outward normal of an implicit surface function is defined as -> V4> N = —— and it. is defined over the entire domain except where \V(f>\ — 0, or where 4> I S not differentiate. This thesis uses finite difference techniques on a cartesian grid to calculate derivatives, such as those in V<f>. 3 Chapter 2 Fast Marching Methods 2.1 Background 2.1.1 Fast Marching Method The Fast Marching Method (FMM) was originally described in [16] and inde-pendently described in [11]. It is a numerical scheme which approximates the solution to the non-linear Eikonal equation [12]. This equation is of the form ||Vu(x)|| =c(x) inQ, c{x)>0 u = b(x) on dfl where fl represents the domain. The boundary is represented by dfl. In the case of the FMM the boundary is the interface and it is fixed throughout the entire calculation. The cost function is c(x) and b(x) is the boundary value function and they are both the given data. The cost function is defined over the whole domain, while the boundary function is only defined on dfl [12]. Equation (2.1) is defined in a continuous domain. The fast marching method uses a discrete method to solve a continuous problem. This discrete method is outlined in Algorithm 1 and it is called Dijkstra's algorithm [3]. In the algorithm the continuous domain fl is discretized into a graph, Q, consisting of nodes, Gn = {xi}, and edges, Qe. The discretized version of the boundary function, dfl, is the set of nodes T C Qn. Each node in Qn has a set of neighbors, M(xi) associated with it. Each non-boundary node, X* € Qn \ T , has a cost, C (XJ) , and each boundary node, x$ G T , has a boundary value, 6(XJ). The algorithm first initializes the value, u(xi), of each node to +00. It then sets the value of foreach x^ £ Qn \ T do u(xi) = +00 for each Xj € T do u(x$) = 6(XJ) G<-e„ while Q j= 0 do Xj <— ExtractMin(Q) foreach Xj £ Af(xi) do U(XJ) *— Update (XJ , Af(xj), u, c) Algorithm 1: Generic shortest path dynamic programming algorithm. De-pending on the choice of update function, this algorithm will produce either Dijkstra's algorithm or the fast marching method. Chapter 2. Fast Marching Methods 4 mm ' m% 5 2 (,) (b) (c) t 4 2 3 2 J (fl I It 11 ^mW ^mW 2 V X (8) y Figure 2.1: Dijkstra's method on a undirected discrete graph. each boundary node with the data. A heap more generally, a priority queue, is constructed, Q, which contains all the nodes in the graph. Dijkstra's algorithm removes the node with the minimum value and updates all of its neighbors by using the specified update procedure. The algorithm continues to do this until Q is empty. Thus, a value function is constructed by dynamic programming. It contains the minimum value, u(xi), at each node, a^, This minimum value represents the cheapest path from a boundary node to xt. Two key factors contribute to the effectiveness of the F M M . First, an up-winding scheme is used to provide an approximation of the viscosity solution for the problem. This type of solution is a particular weak solution which is desirable in this case. Second, the scheme can be computed in a very fast man-ner by marching the solution outwards from the boundary. Since an upwinding method is used, the value at each node increases as the solution is propagated outwards from boundary nodes. Dijkstra's algorithm is usually used to solve a single-source shortest-path problem on a undirected graph. In this case, the cost function is defined on the edges, Ge, but in the case of the fast marching method the cost function is defined on the nodes, Qn. Figure 2.1 demonstrates how Dijkstra's algorithm, Algorithm 1, works using its own update procedure as defined in Algorithm 2. In Algorithm 2, e(xi,Xj) represents the edge between nodes Xj and Xj. If the cost function is defined on the nodes then Dijkstra's method on a discrete graph will never update the value of a node. Thus this example defines the cost function on the edges, c(e(xj, Xj)), to demonstrate how the update function can change a value for a node. In Figure 2.1, (a) demonstrates the graph after initialization. Chapter 2. Fast Marching Methods 5 Input: Xj, M(XJ), u, c Output: mmx^(Xi){c{e(xj,xk))+u{xk)) Algorithm 2: Dijkstra's update method with costs defined on the.edges Input: Xj, Af(xj), u, c Output: C(XJ) + min a ; f c_v( a ! 3.) u{xk) Algorithm 3: Dijkstra's update method with.costs defined on the nodes The boundary node is t G T and it has a boundary value, b(t), of 0, and the value, u(xi), of all other nodes is initialized to +oo. The costs to travel between nodes is defined on the edges, (b)-(g) demonstrate the situation after each successive iteration of the while loop. The value of each vertex is displayed inside the node. The black vertices have had their final values assigned and are no longer in the heap Q. The grey vertex in each image currently has the smallest value in the heap and is next to be removed. The grey edges are the current minimum path from the boundary node t to all other nodes. The update of a node is demonstrated from (d) to (e). Node u originally has a value of 9 and its path is t-u. Then node y is popped from the heap and updates its neighbours, so u has a smaller value of 7 and its new path is t-x-y-u. Dijkstra's algorithm is most efficiently implemented by using a min-heap data structure for Q. Each node in the domain is extracted only once from Q, and this means that it takes O(N) operations to traverse all of the N domain nodes. It takes C^logA^) to maintain the heap property because a node is added once, removed once and updated once for each neighbour and each of these operations is O (log AT) complexity in a min heap. Thus, by using a min-heap data structure Dijkstra's algorithm has a complexity of 0(N\ogN). A backpointer array is maintained in order to determine a node's position in the heap in 0(1) time, and this information is needed in order to update the value of nodes in the heap [12]. The chosen update method determines how the lowest value for each node is defined. A one-norm approximation of this value is calculated by using Di -jkstra's update method on a Cartesian grid with four neighbours. The fast marching method is a two-norm approximation of this value. Dijkstra's update method with the costs defined on the nodes is described in Algorithm 3, and only this version of Dijkstra's update method in implemented in the Toolbox. The fast marching update method is described in Algorithm 4. This implemen-tation of the fast marching method update procedure was originally described in [4], but Algorithm 4 is updated so that all dimensions have the same grid size, h, instead of a grid size of 1 which was assumed in [4]. 2.1.2 Signed Distance Function A signed distance function is an implicit surface function with the additional property that | |Vu | | = 1. A n implicit surface function, cj>, can be turned into a signed distance function, u. The signed distance function will have the same zero Chapter 2. Fast Marching Methods 6 Input: Xj, JV(XJ), u, c for A; € 1,.., n do a*; = mm(u(j\fk(xj))), where A/fe(xj) is the set of Xj's two neighbours along the k coordinate a = SortAscending(a) m = n . while C(XJ)2 < Zke{i,..,m}(h(am - ak))2 do m = m — 1 t = ^ ^Sfcafc + ^/m(/ic(xj)) 2 + EfcS;afca; — mEfca2, j , where E indexes are' {1,.., m} Output : £ Algorithm 4: Fast marching update method level set as the implicit surface function but it will have a gradient magnitude of one. The fast marching method can be used to solve for a signed distance function by choosing c(x) = 1 and V = dfl [!}. The set of boundary nodes, T , and their values, b(xi), must also be defined in order to run the fast marching method. In the case of the signed distance function, the boundary function is an interface which consists of the zero level set of the implicit surface function. Thus, the boundary nodes are the nodes which are closest to this interface. These nodes can be identified by checking to see if the node has at least one neighbour with an opposite sign. If the signs are opposite then the nodes are on different sides of the interface and'thus the node is a boundary node. Once a boundary node is identified its distance, d, from the interface is calculated to be used as its initial value. A n example of a boundary node with two neighbours on the other side of the interface can be seen in Figure 2.2. There are several other cases where the neighbours are in other positions; to see these examples consult [1]. In order to calculate the distance from a boundary node, XQ, to the interface, first a distance, Si, must be calculated for each neighbour Xi of xo which is on the other side of the interface. The node xo has up to 2n neighbours, x\, • • • ,xin, where n is the number of dimensions of the domain. The distance, Sj, is defined as the distance from XQ to the intersection point of the interface, x^, on the line connecting the two grid points. The value for Si is approximated by comparing slopes of </> along the line connecting the two grid points. The slope of (f> for the line segment from XQ to its neighbour x^ is compared to the slope of <j> for the line segment from XQ to the intersection point, ii in order to solve for Sj. Here is the equation comparing the two slopes 4>(xj) - <J>(XQ) _ 0 - <j>(x0) Xi - x0 (Si + x 0 ) - x 0 and here is the same equation rearranged to solve for s» _ h(j>(xo) <f>(x0) - <t>{Xi) Chapter 2. Fast Marching Methods 7 interface^^^^ / , Si = Si / \ ^ / ^ S2 = S2 x2 xx %2 Figure 2.2: Two neighbors of the boundary node xo are on the other side of the interface. where <fi represents the implicit surface function and h represents the grid spac-ing. After all Sj values are calculated the minimum value in each dimension 5j is kept since only the distance to the closest point on the interface is needed. Here is the equation which approximates the gradient of the distance function « i ) + G2) + + Gn) 1 where the right side is 1 and the left side represents an upwind approximation to the gradient. Here is the equation rearranged to solve for d ^ _ sis2 ••• sn y/(s2 • • • sn)2 + (sis3 • • • sn)2 H h (si • • • Sn-1)2 repeat the equation for all boundary nodes and then run the F M M . The fast marching method returns a function u which represents a distance function but not a signed distance function. Fortunately, the sign at each node is already known because the interface moves less than a grid size, h, throughout the calculation so the nodes in the signed distance function have the same sign as those in the initial implicit surface function, <f>. Thus, each node in the distance function must be multiplied by its sign in 4> in order to calculated its signed distance value. Chapter 2. Fast Marching Methods 8 2.1.3 Extension Velocity One of the most fundamental level set methods solves an initial value P D E for motion in the normal direction. The equation for motion in the normal direction is Mt,x)+F(x)\\V4>(t,x)\\=Q (f>(x, t = 0) given It is possible that the original speed function, F, is only defined on the interface and needs to be extended onto all level sets. It is also possible that the speed function does not move the level set representation in a nice manner. The effect of the speed function could cause the representation to develop a large or small gradient magnitude, so <f> becomes very steep, flat or non-monotone.- In both cases the level set equation can be restated as 4>t{t,x)+Fext{x)\\V4>{t,x)\\=0 Fext = F on x) = 0 where Fext is the extension velocity. Extension velocity is a common technique used to extend the definition a speed function from the zero level set onto neighboring level sets [1]. This speed function is a parameter for a level set method calculation. The only restriction on the extension velocity is that the value of the extension velocity should be equal to the value of the original speed function at the zero level set. This restriction is set to ensure that the motion of the implicit surface is the same under (2.2) or (2.3). The stipulation put on the extension velocity in this implementation is VFext-Vu = 0 2 Fext = F on <p(t, x) = 0 where Fext is the extension velocity and u is a signed distance function of the implicit surface function, <j>. The speed function, Fext, can not change in the normal direction of the interface but can only change tangential to the inter-face. Thus, Fext is constant along the normal direction of the interface, so <fr's isosurfaces will maintain their original spacing while being evolved by (2.3). In this implementation, the fast marching method is used to construct an extension velocity and a signed distance function for an initial interface at the same time [1]. The signed distance function is calculated as was described in the previous section, but an additional function is added to the update method in order to calculate the extension velocity. Firstly, extension velocity values, Fext, must be calculated for boundary nodes. A method for identifying these boundary nodes was described in the previous section, as well as a method for calculating the distance, d, between the boundary node and the interface. There are two ways to calculate the extension Chapter 2. Fast Marching Methods 9 velocity value at boundary nodes. The first way returns the original speed, F, defined at the update node. This method must be used if the original speed is only defined on grid nodes. The second method assumes that a function for the original speed is defined throughout the entire domain or at least continuously near the interface. In this case, a weighted average of speed values is taken for the points Xi which were used to compute d. Here is the extension velocity equation which is equivalent to solving (2.4) locally on a grid cell Q = ,'Fext(xQ) - -F(xi) Fext(x0) - F(xn)\ ( d Si sn J \Si t where Fext(xo) is the extension speed for the point XQ and &i is the point where the interface intersects with the Cartesian grid lines as illustrated in Figure 2.2. Here is the equation when solved for Fext(xo) j[F{x1) + j*F{x2) + --- + %F(xn) Fext{Xo) = 1 , 1 , _ 1 ~T + 5T H + ~T Thus, the boundary conditions for (2.4) are defined and next the method for updating extension values for nodes not on the boundary needs to be defined. The extension values are updated at the same time as the distances are updated in the fast marching method. These extension values are calculated by using the same nodes as are used in the update of (f>. These update nodes, Xi where i € {1, • • • ,n}, are defined in Algorithm 4 as the nodes at which the a, values are chosen when XQ is being updated. A n upwinded gradient for the signed distance value and the extension velocity value is chosen at each node to solve (2.4). Here is the upwind version of (2.4) where gradients are taken assuming that there is a node in every dimension which is used to update (p u(xi) - u(x0) u(xn) - u(xp) h h (xi) - Fext(x0) Fext(xn) - Fext(x0) h h = 0 and here is the equation rearranged to solve for .Fext^o), the extension value for the point XQ, p , s Fextjx^jujx!) - u(x0)) H h Fext(xn)(u(xn) - u(x0)) e x t [ X 0 ) ~ (u(Xl) - u(x0)) + ••• + (u(xn) - u(x0)) 2.2 I m p l e m e n t a t i o n This section discusses some general implementation choices. It also states and describes how to use the three methods which were added to Ian Mitchell's Level Set Toolbox [6]. Chapter 2. Fast Marching Methods 10 2.2.1 Implementation Choices The purpose of this implementation is to add new functionality to Ian Mitchell's Level Set Toolbox [6]. Previously, the toolbox only had functions to solve time-dependent Hamilton-Jacobi PDEs, but this new functionality allows some static Hamilton-Jacobi PDEs to be solved. Since the toolbox is implemented fully in Matlab, the new fast marching method was also implemented in Matlab so that the portability of the toolbox would not be lost. Unfortunately, this implementation was extremely slow and even optimizing the code with help from Matlab's profiler could not speed up the implementation very much. Thus, the code had to be re-implemented in C and called from Matlab through a M E X interface. This greatly increased the speed of execution. A n option was to only implement the heap in C, because in Matlab the heap was the slowest part of the code. This option was not chosen because Dijkstra's method is heavily dependent on heap functions. In Algorithm 1 each ExtractMin call requires one heap function call and each Update call requires up to four heap function calls. Thus, 5iV heap function calls would be made for a graph with N grid nodes, so making these calls through M E X interfaces would have added too much overhead to program execution. 2.2.2 Dijkstra This method implements Algorithm 1 and it allows for a specific update function to be chosen. value = d i j k s t r a ( g r i d , bdryData, bdryMask, cost , update_func) The g r i d parameter is a grid structure as defined in Ian Mitchell's Level Set Toolbox [6]. Each dimension in the grid must have the same gr id .dx value. The bdryData parameter is the array defining costs at boundary nodes and the bdryMask is a boolean array defining the boundary nodes. The cost parameter is an array which defines the cost to travel to a node from one of its neighbors. A l l three matrices must have the same size as defined in the grid structure. Finally the update_f unc string specifies the desired update function to be used. The fast marching method, Algorithm 4, can be specified by the string 'fmm', and Dijkstra's update procedure, Algorithm 3, can be specified by the string ' d i j k s t r a ' . 2.2.3 Signed Distance This method computes a signed distance function for a given implicit surface function. The method first determines which nodes are boundary nodes and calculates the value at these nodes. The signedDistan.ee function then calls the d i j k s t r a function and passes in the boundary nodes and boundary values as parameters. The d i j k s t r a function computes a signed distance representation of the implicit surface function. The specifics of the method are explained in Section 2.1.2. Chapter 2. Fast Marching Methods 11 signed_dist = signedDistance(grid, in i tFunc) The g r i d parameter has the same constraints as in d i j k s t r a . The in i tFunc parameter defines an implicit surface function. This array must have the same dimensions as defined in the grid structure. 2.2.4 Extension Velocity and Signed Distance This method uses a fast marching method to construct a signed distance function and an extension velocity for an implicit surface function and an initial speed function. The specifics of the method are discussed in Section 2.1.3. [ex t_veloc i ty , signed_dist] = ex tendVeloc i ty(gr id , i n i tFunc , speedOption, speed_func) The g r i d parameter has the same constraints as in d i j k s t r a . The in i tFunc parameter defines an implicit surface function. This array must have the same dimensions as defined in the grid structure. The speedOption parameter speci-fies whether the speed at boundary nodes will be calculated at only those nodes (1) or whether they will be averaged by speeds from the surrounding interface (2) . The difference between the two choices is discussed in Section 2.1.3. F i -nally, the speed_f uxic parameter must be defined. This parameter can either be a function handle or an array. If the parameter is a function handle then it refers to a Matlab method which will calculate the speed at a domain node. This func-tion handle has one parameter, gr idPoint , a vector representing a point in the domain. The prototype for speed_func is speed = speed_func(gridPoint). If the parameter is an array then it defines the speed at each grid node and the array must have the same dimensions as defined in the grid structure. In the case when speedOption is 1 then both types of speed_func parameters can be passed in, otherwise only a function handle can be used as a parameter for speed_func. 2.3 Examples The purpose of this section is to use some example problems to test the imple-mentation of the methods described in Section 2.2. 2.3.1 Shape From Shading This section describes the function fmmPublish/examples/shapeFromShading. In general, the shape from shading method allows the shape of a surface to be determined from a shaded image of the surface. In this particular case, an Eikonal equation is being solved to determine a viscosity solution for a Lam-bertian surface illuminated by a single distant light source. The test problems Chapter 2. Fast Marching Methods 12 0 0 D 0 . E t M 0 . B 0 J 1 Figure 2.3: Surface (left) and contour (right) plots for the smooth solu-tion of the shape from shading problem. They are generated by the call shapeFromSliading(41, 1, 'fmm')- The domain of this example is 41 by 41 nodes. The fast marching method is used as the update method to solve the problem. were given in [8], while their analytical solutions were given in [9]. The Eikonal equation is defined on a two dimensional domain from [0, +1]2. The cost function is c (x, y) = 27ry [cos (27rx) sin (2wy)]2 + [sin (27rx) COS (2iry)]2. Boundary nodes are defined on the exterior of the domain as well as on five interior nodes. The interior boundary nodes are defined at the following points in the domain: (§, | | , (f, § } , ($, f ) , (§ , f ) , (|, | ) . The value of the exterior boundary nodes is always zero. The boundary value at the interior nodes depends on the type of solution desired. A smooth solution defines the values as [+1, +1, —1,-1,0] respectively, with the resulting analytical solution T (x, y) = sin (27rx) sin (27ry). The continuous but not smooth solution defines the values as [+1, +1, +1, +1, +2] respectively, with the resulting analytical solution T{x,y)={ max | sin(2-7rx) sin(27ry) | \ 1 + cos(27rx) cos(27ry) J [ | sin(27rx) sin(27ry) |, if \x + y- 1| < F-2 /1 < 3, otherwise The smooth solution for the shape from shading problem viewed in Figure 2.3 is generated by the function call below. This function call also generates a continuous but not smooth solution and its results can be viewed in Figure 2.4. The domain for both figures is 41 by 41 nodes and both examples use 'fmm' update procedure. They both demonstrate the use of d i j k s t r a . Chapter 2. Fast Marching Methods 13 Figure 2.4: Surface (left) and contour (right) plots for the continuous but not smooth solution of the shape from shading problem. They are generated by the call shapeFromSliading(41, 2, 'fmm'). The domain of this example is 41 by 41 nodes. The fast marching method is used as the update method to solve the problem. [value, g] = shapeFromShading(numNodes, isSmoothSolution, update_func) : The function above demonstrates the fast marching method by using it to solve the shape from shading problem. The numNodes parameter allows the user to specify the number of grid nodes in each dimension. This parameter has a pre-condition of numNodes = 4m — 1, where m > 1, which is applied to ensure that the interior boundary conditions occur at nodes. The isSmoothSolution pa-rameter specifies whether the smooth solution is desired (true) or whether the merely continuous one is desired (false). Finally the update.func parameter specifies the desired update function to be used in method. The fast march-ing method can be specified by 'fmm' and Dijkstra's update procedure can be specified by ' d i j k s t r a ' . 2.3.2 Extension Velocity of a Circle This section describes the function fmmPublish/examples/circleExtensionVelocity. The purpose of this function is to demonstrate how an extension velocity and a signed distance function can be calculated from an initial function. The initial function chosen in this example is the signed distance from a circle multiplied by a factor of 5. The circle has a radius of 1 and is centered at [0,0]. The initial function is specified by the following equation 4>(x,0) = 5( | | a - | | 2 - l ) . (2.5) Thus, the initial function will be a cone with steep walls and a gradient mag-nitude of 5 as demonstrated in Figure 2.6. The domain of this initial function is Chapter 2. Fast Marching Methods 14 Figure 2.5: Original (left) and extension (right) velocity func-tions of Equation (2.5). The images are generated by the call c i r c l eEx tens ionVeloc i ty (41 , 2). Figure 2.6: Initial (left) and signed distance (right) functions for the im-plicit surface function in (2.5). The images are generated by the call c i r c l eEx tens ionVeloc i ty (41 , 2). Note the vertical scale. calculated in a 2-dimensional space of [—1.5,+1.5]. The original speed at each grid point Xj is the two-norm of that position, F(x;) = ||x»||. A n extension velocity is not allowed to change normal to the interface and the original speed doesn't change tangential to the interface. The original speed is denned as 1 on the entire interface which is a circle with a radius of 1, so the extension velocity will be a flat plane at 1. The signed distance will be a shallower cone than the initial function and it will have a gradient magnitude of 1. The extension velocity can be seen in Figure 2.5 and the signed distance for the problem described above can be seen in Figure 2.6. They are gener-ated by the c i r c l eEx tens ionVeloc i ty function which demonstrates the use of extendVelocity. [ex t_veloc i ty , s igned_dist , gr id] = c i rc leExtens ionVeloci ty(nodes , speedOption) : The numNodes parameter allows the user to specify the number of grid nodes in each dimension. The speedOption parameter specifies whether the speed at Chapter 2. Fast Marching Methods 15 boundary nodes will be calculated at only those nodes (1) or whether they will be averaged from speeds from the surrounding interface (2). See Section 2.1.3 for further details about the speed options. 2.3.3 Movement of a Circle in Normal Direction This section describes the function fmmPublish/examples/normalCircle. The purpose of this function is to demonstrate how to solve a level set equa-tion, for example a time dependent Hamilton-Jacobi equation, while updating its velocity with an extension velocity. The Level Set Toolbox is used to handle the time dependent P D E . This section demonstrates the capability of combining the code from this thesis with the features already in the Toolbox. In this example, an initial function, <j>, is five times a signed distance to' a circle with radius 1 centered at [0,0]: (p(x, 0) = 5(||x||2 — 1). This initial function will be a cone and its walls will have a gradient magnitude of 5. The domain of this initial function is calculated in a 2-dimensional space of [—1.5,+1.5]2. The level set equation, (2.2), moves the initial function in its normal direction. Thus, the surface at the zero level set will be a circle and will expand into a larger circle which is still centered at the origin. In this case, the initial speed function at each grid point is the two-norm of that position, F(x) = \\x\\2- Thus, without the extension velocity the initial cone will move faster the further away it is from the origin. This will result in the cone becoming more shallow. The result with the extension velocity will be slightly different. A n extension velocity is not allowed to change normal to the interface and the original'velocity doesn't change tangential to the interface. The original velocity is defined as 1 on the entire interface, so the extension velocity will have a value of 1 throughout the entire domain. Thus, the walls of the cone will move at approximately the same rate and the gradient magnitude of the walls should remain near 5. The bottom of the cone will flatten out. The different results which the extension velocity causes in the normalCircle example can be viewed in Figures 2.7 and 2.8. Both figures demonstrate that without an extension velocity the gradient magnitude decreases with time and with an extension velocity the gradient magnitude remains the same for most nodes. [g, speedStore] = normalCircle(useExtVel, accuracy, numNodes) The method above computes an implicit surface function at four timesteps and at each timestep the contours of the function are plotted. The implicit functions for each timestep are returned in a cell vector called speedStore. The grid which represents the domain of the implicit function is also returned. This grid parameter is a grid structure as defined in Ian Mitchell's Level Set Toolbox [6]. useExtVel is a boolean parameter which is true if an extension velocity will be used and false otherwise. The accuracy parameter specifies the degree of accuracy used to solve the level set method. The choices are 'low', 'medium', Chapter 2. Fast Marching Methods 16 I - 0.125 1.0 1.0.125 Figure 2.7: Contour plots for a circle moving in its normal direction. The left im-age demonstrates the level set method run with its regular speed, F(x) = ||a-||2. It is generated by the call normalCirc le (0 , 'medium', 41). The contours become further apart as the cone shaped surface becomes shallower. The right image demonstrates the level set method with an extension velocity. It is gen-erated by the call no rma lCi rc l e (1 , 'medium', 41). The contours stay the same width apart because the whole surface is moving in its normal direction at the same rate. ' h i g h ' and 'veryHigh ' . The numNodes parameter allows the user to specify the number of grid nodes in each dimension. Chapter 2. Fast Marching Methods 17 —I-0.126 - " i - o a • I - 0.37B 1000 1200 1400 1B00 Nona Numbar 1200 1400 1600 Figure 2.8: Ordered gradient magnitude plots for the surface at different times. In the plots the nodes are sorted by gradient magnitude. The left image demon-strates the level set method run with its regular velocity. It is generated by the call normalCirc le(0 , 'medium', 41). The most common gradient mag-nitude decreases as the surface gets shallower. The right image demonstrates the level set method with an extension velocity. It is generated by the call normalCirc le(1 , 'medium', 41). The most common gradient magnitudes are all 5. 18 Chapter 3 Reachable Manifold Sets on a 3.1 Background 3.1.1 Differential Algebraic Equations A D A E is a differential equation which is comprised of algebraic (no derivative is present) and differential (derivatives are present) terms [2]. Here is the general implicit form for a D A E : where dF/dy' may be singular or nonsingular depending on the index of the D A E . The index of a D A E refers to the number of times differentiation must be carried out on the system until an ordinary differential equation is obtained for y. In other words, the number of differentiations required to solve uniquely for y' in terms of y and t [2]. This research only deals with DAEs of index-1, since we need to differentiate the algebraic term in order to specify the velocity field. More specifically this research deals with semi-explicit DAEs of index-1 or in other words an O D E with constraints where the constraints represent a manifold on which the system state evolves. The general form for this type of where x represents the differential variables and z represents the algebraic vari-ables. In this form it is easy to identify g as the constraint function since its equation does not contain derivatives. 3.1.2 Closest Point Method The closest point method [10] allows a P D E to be extended off the manifold onto the entire domain. The problem, which is originally only defined on the surface, is embedded onto the whole cartesian grid so that standard level set techniques can be used to solve it. The solution on the manifold can then be obtained by only considering the results along the manifold. This approach is F{t,y,y')=0 D A E is x' — f(t,x,z) 0 = g(t, x, z) (3.1) Chapter 3. Reachable Sets on a Manifold 19 an alternative to either parameterizing the constraint manifold or triangulating the constraint manifold in order to solve the P D E . The embedding approach is accurate and the P D E is easy to solve numerically once the embedding is complete because there are a lot of well known techniques for solving PDEs on Cartesian grids. Fortunately, the embedding approach used in the closest point method is easy to implement as well. The embedding approach creates an extended function, which involves re-placing surface gradients, defined only on the manifold surface, with standard gradients which are defined on the entire computation domain, A require-ment for this procedure is that the embedded P D E must be a natural extension of the surface P D E . The simplest way to meet this requirement is for the gradi-ent to be constant along the direction normal to the surface. This means that the gradient will only change in the direction of the surface and no artificial rates of change will be added to the P D E . The closest point method uses an extended function which is created by composing the original surface function onto a closest point operator. A simple and accurate way to define this closest point operator, CP(x), for a point, x, is to have the operator return the closest point to x on the surface. This operator is used before each timestep in the level set solve to reinitialize the implicit surface function, <j>, as 4>(CP(x)). This will ensure that for a point on the surface, XQ, all points, xn, along its normal direction will have the same function value, 4>{XQ) = 4>{xn). Thus, the gradient of 4>(CP(x)) will be constant along the normal direction and will only be able to change in directions tangential with the surface. This closest point operator will create a well defined and smooth function near the surface, but may have discontinuities away from the surface where a point is equidistant from two points on the surface. The greater the local curvature in the manifold the closer the discontinuity will be to the manifold, but the discontinuity should be at least several grid cells away from the manifold in order to achieve accurate results. The closest point method allows for a variety of difficult surfaces to be represented such as open and closed surfaces in any codimension. The codimension of a surface with dimension p in a domain of dimension n is n — p. We only study codimension 1 surfaces in this thesis, but the closest point method should work for higher codimensions. Research in [15] describes a method for calculating the closest point operator. This method uses Voronoi Regions to calculate the operator. Analytic calcu-lation, using standard numerical optimization techniques or using a tree-based algorithm for triangulated surfaces are options suggested in [10] in order to cal-culate the closest point operator. The research in the thesis solves a constrained non-linear optimization problem in order to determine the closest points. 3.1.3 Reachable Sets Model checking is a method used to verify or validate complex engineering sys-tems. A common form of model checking is to verify whether a model in evo-lution can enter into an unsafe state. This form of model checking involves the Chapter 3. Reachable Sets on a Manifold 20 computation of a reachable set of which there are two types: forwards reachable sets and backwards reachable sets. Computing forwards reachable sets involves evolving an initial set of states forward in time in order to determine the set of states reached by their trajectories. Computing backwards reachable sets involves declaring a set of target states and then determining the set of start states whose trajectories reach that target set. In [14], a method was devel-oped to solve backwards reachable sets using level set methods. This method expresses the reachable set problem as a Hamilton Jacobi equation. Level set methods are a well known technique for solving the Hamilton Jacobi Equation (HJE), thus they can be used to solve backwards reachable sets. In order to solve for the backwards reachable set, a continuous system with dynamics, x = f(x), is defined in order to move a set of initial conditions: The backwards reachable set is the subset of these initial conditions that are driven into the target set, TQ. The state of the system is defined by x and the target set is defined as a zero sublevel set To = {x € SR" | <j>o{x) < 0} of a scalar function (fi(x,0) = <fio{x)- The backwards reachable set is defined as T (T ) over a finite horizon r < oo, the trajectory of the system is defined by a;(-) and the state of the trajectory at time r is defined by X{T). Thus, a formal definition of the reachable set, T ( T ) , is the set of a;(0) such that X(T) € %. It can also be defined in implicit surface notation as a zero sublevel set T( r ) = { i £ S n | <t>(x, -T) < 0} where 4>(x,t) is the solution to the following H J E which is solved backwards in time from t = 0 to t = —r < 0 0 = Dt(t>(x,t) + f(x,t)-V(j>(x,t) (3.2) The reachable tube is another object sought from this algorithm. The reach set is the set of states occupied by trajectories at a specific point in time r; alternatively, the reach tube is the set of states that have been traversed by the trajectories from t = 0 to t = T [5]. More formally, it is the set of x(0) such that x(s) € % for some s € [0, r]. The following H J E computes the reach tube 0 = A<A(x, t) + min[0, f(x, t) • V<j>(x, t)\ (3.3) The minimization with zero in the second term constrains the temporal derivative to a positive sign. This restriction keeps the reachable tube from shrinking as t —» — oo. This research calculates reachable tubes for the two-dimensional toy problem in Section 3.3.1 and it calculates reachable sets for the three-dimensional toy problem in Section 3.3.2. Chapter 3. Reachable Sets on a Manifold 21 3.2 Implementation 3.2.1 Level Set Methods The implementation is accomplished by using Ian Mitchell's Level Set Toolbox for Matlab [6]. The toolbox uses upwinding schemes to solve H J E which move under a convective field such as (3.2) and (3.3). The solutions displayed in this work will use second order accuracy unless otherwise stated. The convec-tive field is approximated by termConvection, the upwinding spatial derivative is approximated by upwindFirstEN02 and the temporal derivative is approx-imated by odeCFL2. The function termRestrictUpdate is used in order to restrict the temporal derivative to a positive sign such as required in (3.3). If the restriction is used then a reach tube is calculated, otherwise a reach set is calculated. 3.2.2 Differential Algebraic Equations The D A E specifies the velocity field required by termConvection. Since the D A E required for this research is an O D E with constraints, the constraint term must be differentiated with respect to t in order to specify a velocity field for all variables. Here is the velocity field, which is derived from differentiating (3.1) and rearranging to solve for z', f(t,x,z) _ dg(t,x,z) \ _ dg(t,x,z) dx J \ i > / dt dg(t,x,z) dz 3.2.3 Closest Point i This implementation of the closest point operator assumes that the constraint surface is represented by an implicit surface function. That implicit surface function is then used to compute a signed distance function using the methods in Chapter 2.1.2. Algorithm 5 demonstrates how the closest point operator is computed us-ing the signed distance function and the implicit surface function. Table 3.1 describes the parameters used in the closest point operator algorithm, while Tables 3.2 and 3.3 describe the inputs and outputs for the functions used in the closest point operator algorithm. The algorithm takes in a signed distance representation of the manifold, signedDist. It then sorts, SortAscending, the nodes in this representation depending on the absolute value of their distance.1 The nodes are then traversed 1 A topological sort, where each node occurred after its neighbour of lowest value, would be faster than a strict sort by distance but in practice the optimization routine dominates running time. Thus, the fraction of time saved by using a topological sort is negligible. dx H dz ~di Chapter 3. Reachable Sets on a Manifold 22 Input: signedDist, constraint sorted = SortAscending(Abs(signedDist)) foreach Xj £ sorted do if IslnterfaceNeigh(xj, signedDist) then x — TravelGradient(xj, signedDist) else xs = SmallestNeigh(a-j, signedDist) x = closestPoint(x s) x<?p = CpOptimizedj, x, constraint) closestPoint(XJ) — Xcp Output : closestPoint Algorithm 5: Algorithm which computes the closest point operator. Parameter Name Type Description signedDist Array Signed distance representation of the manifold constraint Function handle Functional representation of the implicit surface function for the manifold sorted Vector Sorted list of nodes depending on their distance to the manifold closestPoint Cell vector with ar- Closest point on the manifold for ray elements each grid node Table 3.1: Description of parameters used in the closest point operator algo-rithm. in their sorted order from the smallest distance to the manifold. If the node is adjacent to the interface, Islnterf aceNeigh, then the algorithm takes steps, on the order of Ax, along the gradient until the interface is reached. This step was needed in order to give an accurate initial point x for the optimization problem. To compute the initial guess for nodes which are not adjacent to the interface first the neighbor which is closest to the interface, SmallestNeigh, is found. The initial guess becomes this neighbor's closest point, closestPoint, which has already been computed. The initial guess, x, is used in a nonlinear constrained optimization problem, CpOptimize, which finds the.closest point on the manifold. Here is the optimization problem: min | | x c p - X i | | 2 subject to 0 =constraint(xcp) where xcp is the closest point on the interface for node Xj and constraint (x) is the function representing the manifold. The constrained optimization problem is solved using Matlab's fmincon function. The closest point function is composed with </> to reinitialize ^ as an em-Chapter 3. Reachable Sets on a Manifold 23 Function Name Inputs SortAscending A List of real numbers Is lnterfaceNeigh A grid node, signedDist TravelGradient A grid node, signedDist SmallestNeigh A grid node, signedDist CpOptimize A grid node, an initial guess for position of closest point, constraint Table 3.2: Description of functions and their input parameters used in the closest point operator algorithm. Function Name Outputs SortAscending List sorted in ascending order Is lnterfaceNeigh True if the node is beside the interface, false other-wise TravelGradient Estimated position of closest point on manifold along the gradient SmallestNeigh The neighboring node with the smallest signedDist value CpOptimize Closest point on the manifold to the grid node Table 3.3: Description of functions and their output parameters used in the closest point operator algorithm. Chapter 3. Reachable Sets on a Manifold 24 bedded surface. The closest point function returns a point, xo = CP(xn), on the manifold, but this point is not necessarily a grid node. Thus, the function value of <p(xo) must be interpolated in order to reinitialize <f> at each node. The method used for this interpolation calculation is Matlab's interpn with cubic accuracy. 3.3 Examples The validity of the method developed in this research is demonstrated in some simple toy problems. Matlab's D A E solver can be used to generate sample trajectories for semi-explicit DAEs. of index-1. A trajectory's position at specific times can be calculated from these trajectories and compared to the implicit surface calculated by the level set implementation. These problems are easy to conceptualize and visualize thus the results of our level set implementation can be clearly grasped. A l l examples are calculated on a Linux box with a Trison Pentium 4 3GHZ processor, I G R A M and 80G Harddrive. The software used was Suse linux 10.1, Kernel version 2.6.16.27-0.9-smp, Matlab 7.3.0.298 (R2006b) and G C C version 4.1.0. 3.3.1 Two-dimensional D A E Toy Problem The D A E for the two-dimensional toy problem is as follows xi = sm(nx2) ^ " X2 = 1.5 + Xi where the first equation is the constraint. In order to derive the velocity field the constraint is differentiated. The velocity field is as follows ii = vi = ivy? cos(nx2) X2 = v2 = 1.5 + Xi • The constraint surface is a sine curve in the vertical direction, while the dynam-ics move along the sine curve faster on the right side of the curve than on the left. The implicit surface representing the initial reach set is a hyperplane with a normal of [0,1] which passes through the origin. The problem was defined on a domain of [—1.5,-1-1.5] x [—0.2, +3.7], while the distance between grid nodes is defined as 0.02. The P D E is executed from 0 to 3 time units and reachable sets are calculated every 0.25 time units. The toy example is calculated using the closest point method and without using the method in order to demon-strate the purpose of the method. It takes'620 seconds to solve the problem using the closest point method and 130 seconds to solve it without using the Chapter 3. Reachable Sets on a Manifold 25 closest point method. The closest point method takes additional time for calcu-lation because it has to interpolate the closest point values for all nodes at each timestep. Figure 3.1 demonstrates that the implicit surface which uses the clos-est point method is perpendicular to the constraint surface, while the implicit surface which does not use the method is not perpendicular to the constraint. The solution using the closest point method is calculating the reachable set on the manifold, while the solution without the closest point method is calculating the reachable set for the whole domain. The calculation solely on the mani-fold makes it much clearer to view where the reachable set intersects with the manifold and it makes it more numerically stable as well. The closest point method also improves the quantitative solution as time pro-gresses in the calculation. This improvement is shown by using Matlab's D A E solver on (3.4) to very accurately estimate the position of the trajectory at each timestep, then the implicit surface value is interpolated at each position of the estimated trajectory. The values at each timestep can be seen in Table 3.4, and since the reach set interface is defined as the zero level set all values should be zero. In the first couple of timesteps the method which does not use the closest point method has a slight advantage, but midway through the calculation that method loses its advantage and becomes far more inaccurate. The closest point method allows the level set calculation to maintain better accuracy throughout the calculation. A convergence plot for t = 3.0 time units, in Figure 3.2, shows superlinear convergence for the two dimensional toy example. The numerical values for the rates of convergence are calculated in Table 3.5. We are not sure why the convergence rate becomes negative when 600 grid nodes are used. These values average out to a convergence rate of 1.3910, which is not quite the quadratic rate expected from the second order accurate scheme used in this research. This may be due to a slight mis-approximation in the closest point calculation where some closest points are slightly off the constraint surface. 3.3.2 Three-dimensional D A E Toy Problem with Codimension 1 The D A E for the three-dimensional toy problem is as follows X \ = sin(7ra;2). x 2 = x 3 X3 = - X 2 where the first equation is the constraint and the last two equations are the differential equations. In order to derive the velocity field the constraint is differentiated. The velocity field is as follows Chapter 3. Reachable Sets on a Manifold 26 Timestep (time units) No Closest Point Closest Point No C P - C P 0.00 0.0000 0.0000 0.0000 0.25 0.0000 0.0009 -0.0009 0.50 0.0001 0.0008 -0.0008 0.75 0.0005 0.0010 -0.0005 1.00 0.0010 0.0017 -0.0007 1.25 0.0069 0.0129 -0.0060 1.50 0.0509 . 0.0109 0.0400 1.75 0.1723 0.0157 0.1566 2.00 0.2439 0.0076 0.2363 2.25 0.2338 0.0098 0.2240 2.50 0.2565 0.0116 0.2450 2.75 0.2571 0.0125 0.2446 3.00 0.2634 0.0217 0.2418 Table 3.4: Interpolated implicit surface value for each timestep of the two-dimensional toy D A E problem. Values using the closest point method and values without using the method are compared. Number Grid Nodes Error Rate 75 0.04923 -150 0.0078 2.6580 300 0.002721 1.5193 600 0.002729 -0.0042 Table 3.5: Convergence of error for the last timestep, t = 3.0 time units, for the two dimensional toy example. -Chapter 3. Reachable Sets on a Manifold 27 ii = v\ = 7ri>2 cos(7ra2) X2 = V2 = X3 X3 = V3 = - X 2 The constraint surface is a sine curve in the X2 direction, while the dynamics move along the sine curve in a clockwise circle in the X2 versus £ 3 plane. The initial reachable set is a cuboid (a rectangular box) which has an in-finite length in the x\ dimension. The upper corner of the rectangle is at (xi, 0.25,0.75) and the lower corner of the rectangle is at (x\, —0.25,0.25). The problem was defined on a domain of [—1.0,+1.0]3, while the distance between grid nodes is defined as 0.05. The P D E is executed from 0 to 27T time units and reachable sets are calculated every 7r /4 time units. The toy example is cal-culated using the closest point method and without using the method in order to demonstrate the purpose of the method. It takes 848 seconds to solve the problem using the closest point method and 392 seconds to solve it without using the closest point method. The closest point method also slightly improves the quantitative solution for the three-dimensional example. This improvement is shown by using Matlab's D A E solver to calculate the position of many sample trajectory at each timestep. These trajectories are initiated from points along the target set, these points have a spacing of 0.02. For each trajectory, the implicit surface value at each position is interpolated. The average values at each timestep can be seen in Table 3.6 and the maximum value at each timestep can be seen in Table 3.7, since the constraint surface is defined as the zero level set all values should be zero. On average at each timestep the solution which uses the closest point method is better than the solution which does not use the method. On the other hand, the maximum value is higher at each timestep when using the closest point method. Chapter 3. Reachable Sets on a Manifold 28 Timestep (time units) No Closest Point Closest Point No C P - C P 0.00 0.0003 0.0004 0.0000 0.7854 0.0103 0.0083 0.0020 1.5708 0.0163 0.0151 0.0012 2.3562 0.0213 0.0191 0.0023 3.1416 0.0256 0.0205 0.0052 3.9270 0.0300 0.0252 0.0048 4.7124 ' 0.0334 0.0316 0.0018 5.4978 0.0371 0.0343 0.0028 6.2832 0.0404 0.0380 0.0024 Table 3.6: The average interpolated implicit surface value for multiple O D E trajectories at each timestep of the three-dimensional toy D A E problem. Val-ues using the closest point method and values without using the method are compared. Timestep (time units) No Closest Point Closest Point No C P - C P 0.00 0.0036 0.0036 0.0000 0.7854 0.0479 0.0483 -0.0003 1.5708 0.0595 0.0645 -0.0050 2.3562 0.0695 0.0757 -0.0062 3.1416 0.0759 0.0799 -0.0039 3.9270 0.0827 0.0963 -0.0136 4.7124 0.0877 0.1014 -0.0137 5.4978 0.0933 0.1158 -0.0225 6.2832 0.0979 0.1273 -0.0294 Table 3.7: The maximum interpolated implicit surface value for multiple O D E trajectories at each timestep of the three-dimensional toy D A E problem. Val-ues using the closest point method and values without using the method are compared. Chapter 3. Reachable Sets on a Manifold 29 (•3 -1.5 -1 -0.5 0 0.5 1 1.5 t = 3 -1.5 -1 -0.5 0 0.5 1 1.5 Figure 3.1: Reachable set plots for the two dimensional toy example. The top image doesn't use the closest point method, while the bottom image does use the method. In both images, the dotted curve represents the manifold, the stars represent the solution calculated using Matlab's D A E solver at intervals of 0.25 time units and the solid curves represent the reachable set solution at intervals of 0.25 time units. Notice how the reachable set is perpendicular to the manifold near the manifold when using the closest point method. Chapter 3. Reachable Sets on a Manifold 30 1 2 S 10 10 10 Grid Size Figure 3.2: Convergence rates for the last timestep, t — 3.0 time units, for the two dimensional toy example. The x-axis defines the number of grid nodes in the xi dimension of the grid, they are as follows 75, 150, 300 and 600. Chapter 3. Reachable Sets on a Manifold 31 t- o 1-1 Figure 3.3: Reachable set plots for the three dimensional toy example in 3D. In both images, the red surface (darker surface) is the manifold, the green curve (lighter curve) is the reachable set on the manifold and the blue dots represent the solution calculated using Matlab's D A E solver. The top image represents the system at the starting time, t = 0 time units, and the bottom image represents the system at the end time, t = 2ir time units. Chapter 3. Reachable Sets on a Manifold 32 t - o -1 0 1 t = 2.3562 -1 0 1 t = 4.7124 -1 0 1 t - 0.7854 -1 0 1 1-3.1416 -1 0 1 t = 5.4978 -1 0 1 t- 1.5708 CD i b i G-| 0 -1 I 1 i .1 -1 0 1 t = 3.927 1 1 1 0 0 b (~) •"el! -1 W -1 -1 0 1 1 - 6.2832 i 1 &\ n i n O U 0 -1 -1 -1 0 1 1 • 0 1 - 0.7654 -1 0 1 t - 2.3562 -1 0 1 t - 4.7124 -1 0 1 t-3.1416 -1 0 1 t- 5.4978 t - 1.5708 CD i 0 0 -1 : : 0 '• '• -1 -1 0 1 t - 3.927 ;• : 1 1 0 0 1 -1 0 1 t- 6.2832 CD - 1 0 1 - 1 0 1 - 1 0 1 Figure 3.4: Plot of the reachable set on the manifold in x2 versus x3 dimensions for the three dimensional toy example. The top image doesn't use the closest point method, while the bottom image does use the method. Chapter 3. Reachable Sets on a Manifold 33 t - o t - 0.7854 t - 1.5708 0. -2 0 2 1-3.1416 CD -2 0 2 t - 5.4978 1 r* -2 0 2 -2 0 2 -2 0 2 t - 0 -2 0 2 t - 0.7854 0 -2 0 2 1-3.1416 -2 0 2 t - 5.4978 -2 0 2 1 - 1.5708 1 r r 0! -2 0 2 t - 3.927 -2 0 2 t - 6.2832 CD -2 0 2 Figure 3.5: Plot of the reachable set on the manifold in s, the horizontal distance along the manifold from the origin in the x2 dimension, versus X3 dimensions for three dimensional toy example. The top image doesn't use the closest point method, while the bottom image does use the method. 34 Chapter 4 Predicting Voltage Instability of a Power System 4.1 Background Voltage instability is an important factor to consider when constructing power systems. In the case of a load variation or an event disturbance the system should be able to maintain an acceptable level of voltage. Failure to maintain this level can result in a brown out, which can damage electrical equipment. A numerical method for predicting voltage instability was introduced in [13], and this method models the voltage dynamics by a nonlinear hybrid automata. The automaton combines continuous voltage dynamics with discrete operations . of the power system. The research in this thesis focuses on numerical methods for solving the continuous component of the voltage dynamics, which can be modeled with a D A E . The following D A E which models these continuous voltage dynamics was introduced in [17] , •, Xj+x'^, , xd-x'dE2+x'Q(E) , „ . T d o E = — E + — E> + E f d ' TEFD = -(EFD - E°FD) - K {EQ(E) - E R } , ( 4 1 } 0 = E'2E2 - {x'Pf - {x'Q(E) + E2}2 = g(E',E;x1), where Ea(E) = ^{x1P)2 + {x1Q(E)+E2}2, . x' = xi+xd, (42) P = P 1 — ± mi Q(E) -Qo + HE + BE2. The physical meaning and values of the variables and parameters are described in Table 4.1. Chapter 4. Predicting Voltage Instability of a Power System 35 Several surfaces are defined for the D A E of (4.1). The first is the constraint surface L = {(£', Efd, E) € SR3 \g(E', E) = 0 } , which consists of a two-dimensional manifold in a three-dimensional space. The second surface Efd,E) e dg_ dE (£ ' ,£) = o } is a singular surface in which the model breaks down because the solution does not hold uniqueness properties. The constraint from (4.2) is differentiated in order to derive the following velocity field TdoE = —E + — + Efd, TEfd = -(Efd - E°fd) - K {EG(E) - Er} , (4.3) E'E'E E • 2{x'Q{E)+E2)-E'2 The denominator for E provides the equation for ^ = 2(x'Q(E) + E2) — E12, and the singular surface occurs where 2(x'Q(E) + E2) — E'2 = 0. In mathematical terms this part of the surface needs to be modeled with a D A E with a index higher than one. 4.2 Closest Point Method used as an Extension Velocity The singular surface S is a challenge: not only is the model inaccurate at S, but as we approach 5 the O D E (4.3) equivalent to the D A E (4.1) becomes i l l conditioned (the term for E blows up) so something must be done to the O D E in order to make it behave in an acceptable manner. Some of the undesirable states of the power system are defined by the unsafe set G = {(E',Efd,E)\E<Ec} The states in G are physically unacceptable for operations because the load bus voltages in those states are too low. Thus, these unsafe states can be used to define the target set, % = G, which is shown in Figure 4.1. The dynamics in the target set can be damped out because they are irrelevant to computing the reach set outside the target set. The following equations are used to damp out the dynamics Chapter 4. Predicting Voltage Instability of a Power System 36 Description Name Value generator voltage behind transient reactance E' field excitation Efd load bus voltage E generator bus voltage EG open-circuit transient time constant T' 5s transmission reactance (two routes) Xi 0.1 d-axis synchronous reactance Xd 1.2 d-axis transient reactance x'd .0.2 time constant of first-order model of A V R T 1.5s nominal field excitation ^fd 1.6 gain constant of first-order model of A V R K 7 set-point value of generator bus voltage Er 1 mechanical input power to generator P 1.0 constant reactive power of load Qo 0.5P m current source of load H 0 impedance load B 0 critical value of load bus voltage Ec 0.7 Table 4.1: Physical Meaning of Variables and Parameters in Nonlinear Hybrid Automaton [13] <t>o < -2e 2e < cj>0 < 0 e < 0 F(x) = H(MF(x) where e = 2h and h represents the grid spacing [7]. This damping function is equal to 1 everywhere outside the target set, so outside the target set the func-tion being damped stays the same. Inside the target set the damping function smoothly and quickly damps out F(x) to essentially zero. The function goes to zero within four grid cells. The implementation in this research uses cj>o(x) not CP(4>o(x)) as input to the damping function. The initial 4>o{x) function is smoother, thus after damping this function it produces better results in the reach tube calculations. Unfortunately, not air of the singular surface is inside the target set, so another technique needs to be used to make the dynamics outside the target set behave in an acceptable manner. The extension velocity method, described in [1] and in Section 2.1.3, extends a speed function, which moves the interface in its normal direction, off an interface. This method maintains the signed distance property for the level sets representing the interface. Several attempts were made to adapt this well known extension velocity method to the case f 0 1 , (00 + e) , 1 . 2 + 2e + 2 7 r S m l e 1 Chapter 4. Predicting Voltage Instability of a Power System 37 of extending velocities off manifolds. This proved difficult and inappropriate for several reasons. . In order to calculate the extension velocity the velocity function had to be converted to a speed function in the normal direction of the implicit surface representing the reach tube. This conversion had to be done at each timestep because the reach tube is evolving at each timestep. After the conversion into a speed function, the extension velocity method was used to extend the speed function two times. It was firstly used to extend the speed function in the normal direction of the constraint manifold and then the resulting speed function was extended in the direction of the reach set implicit surface. This conversion seemed inappropriate for the case of extending off the constraint manifold since the extension velocity was being converted to a speed function in the normal direction of the reach set and not in the normal direction of the constraint manifold. Another reason the conversion was inappropriate is that even though the velocity function and the manifold, do not change over time, a calculation was being done at every timestep to try and extend the velocity function off of the manifold. This calculation at every timestep was an unnecessary amount of overhead. A second attempt was made at extending the velocity field off of the con-straint manifold. This attempt involved using the previously described extension velocity method to extend each component of the initial velocity field separately off of the constraint manifold. This extension only had to be done once before the level set calculation began because the constraint manifold and the initial velocity do not change throughout the level set calculation, and this extension could be done separately on each component of the velocity field because the components are independent of each other, This extension got rid of the sin-gularity problems in the dynamics and maintained the correct dynamics on the manifold. Unfortunately, running these extension velocity functions took a lot of extra time and this extra time seemed wasteful when a closest point function was already being computed to reinitialize the implicit surface function repre-. senting the.reach tube. Thus a new extension velocity method is proposed in this research which makes use of the closest point function. This new method for extending velocities off of the constraint manifold as-sumes that the dynamics are well defined on the manifold but may.have problems somewhere else in the domain. The simple solution proposed in this research is to use the closest point operator,.described in Section 3.1.2, to extend each component of the velocity function off of the manifold. For example, if the dynamics F(x) are unsuitable for the HJE's (3.2) or (3.3) we replace it by . . . Fext(x) = F(CP(x)) •' (4.4) The initial velocity function is parallel to the manifold on the manifold, so off of the manifold the extension velocity function will be parallel to the closest part of the manifold. The dynamics are well behaved on the entire manifold because the singular surface, S, only intersects with the manifold inside the target set and the dynamics inside the target set have already been damped out. Thus, extending these well behaved dynamics off of the manifold will Chapter 4. Predicting Voltage Instability of a Power System 38 create well behaved dynamics throughout the entire domain. Two types of figures are generated to show the validity of the method. Fig-ure 4.2 compares trajectory simulations generated in three different ways by Matlab's D A E / O D E solver. The first method for computing the trajectories uses Matlab's D A E solver with the dynamics described in (4.1). The second method for computing the trajectories uses Matlab's ODE solver with the dy-namics described in (4.3). The final method also uses Matlab's O D E solver but it uses the extension velocity dynamics described in (4.4). These extension velocity dynamics are only defined on a discrete grid, so values off the grid must be interpolated by interpn. In the figure, all three trajectories are practically indistinguishable, and this similarity shows that the extension velocity method does not perturb the dynamics significantly. Figure 4.3 calculates some sample trajectories using the interpolated ex-tension velocity, and the starting points of these trajectories are shifted off of the manifold. These shifted trajectories remain parallel to the solutions on the manifold and this shows that the extension velocity contains dynamics which are constant in the normal direction of the manifold. This consistency in the nor-mal direction ensures that extension velocity is a natural extension of the initial surface velocity. Thus, the extension velocity technique used in this research removes the singularities in the domain, and creates well behaved dynamics throughout the entire domain. 4.3 Implementation The implementation of the power generator problem is similar to the imple-mentation of the two toy examples in Chapter 3. The first difference is that the closest point operator is used to extend the velocities off the manifold as described in the previous section. This extension is only done once during ini-tialization. The second difference is that all derivative approximations are third order accurate. 4.4 Examples The D A E for the Power Generator was described in (4.1) and the O D E which represents the velocity field was described in (4.3). The initial reachable tube is a hyperplane with a normal of (1,0,0) passing through the point (Ec, 0,0). The zero level set of the initial reach tube is a plane spanning the E' and Efd dimensions, and this plane divides the E dimension at Ec. This zero level set can be seen in Figures 4.1 and 4.4. The problem was computed on a domain of [+0.5,+1.7] x [+0.8,+1.8] x [-0.1,+0.7], while the distance between grid nodes is defined as 0.02. The domain values in the Efd dimension were scaled down by a factor of 10 compared to the original problem specifications in [13] and in [17]. The P D E is executed from 0 to 5 time units and reachable tubes are calculated every 0.2 time units. It takes 2806 seconds to Chapter 4. Predicting Voltage Instability of a Power System 39 solve the problem using the closest point method when the closest point function has already been precomputed. The problem cannot be calculated without the closest point method. Figures 4.4, 4.5 and 4.6 demonstrate the reachable tube at every integer time unit of the calculation. These figures show that the closest point reinitialization for the implicit surface function forces the reach tube to remain perpendicular to the manifold throughout the calculation. Figure 4.7 shows several views of the reach tube at the final timestep, t = 5. Figure 4.8 compares the interface for the reachable tube at the final timestep with simulations calculated by Matlab's D A E solver using the dynamics described in (4.1). This figure shows that the level set calculation is doing a poor job calculating the reachable tube. The calculation is overestimating the size of the reachable set, so a lot of safe tra-jectories are being lost. Several ways to improve the quality of the results could be to increase the grid resolution or to increase the accuracy of the calculation. The closest point function could also be slightly inaccurate, so improvements to the closest point method in order to calculate the closest points more accurately could improve the results. Chapter 4. Predicting Voltage Instability of a Power System 40 t-o ::::;:::::::::::::p::::::::::j:::::::::::::p:::::::::-_2 Figure 4.1: Target set and singular set plots for the power generator problem. The top image is a top view of the intersecting surfaces and the bottom image is a side view of the intersecting surfaces. In these images, the red surface is the constraint manifold, the green surface is the interface for the target set, the blue surface represents the target set and the yellow surface represents the singular set. Chapter 4. Predicting Voltage Instability of a Power System 41 / / ..../.. / - i — ~ i f ; : : 1.7 1.B 1.5 1.4 1.3 1.2 0.6 0.8 1 1.2 E 14 1.6 Figure 4.2: Trajectory simulations generated in three different ways by Mat-lab's D A E / O D E solver. The blue trajectories are calculated by Matlab's O D E solver using the dynamics described in (4.3), the red trajectories are calculated by Matlab's D A E solver using the dynamics described in (4.1) and the yellow trajectories were calculated by Matlab's ODE solver using the extension velocity dynamics. The side view (top) demonstrates that the three trajectories are very similar when looking at them head on, but the top view (bottom) demonstrates that the dynamics force the trajectories slightly off of the constraint manifold in the bottom left hand corner of the plot. Chapter 4. Predicting Voltage Instability of a Power System 42 0.7 r o -0.1'— ' 1 1 1 t -0.6 0.8 1 1.2 1.4 1.6 E 0.6 0.8 1 1.2 1.4 1.6 E Figure 4.3: Sample trajectories using the interpolated extension velocity. The top image is a side view of the trajectories on the manifold and the bottom image is a top view of the trajectories. The middle set of trajectories is on the manifold and the outer two sets of trajectories are shifted off of the manifold by 0.1 units in the E' dimension. The trajectories remain parallel to the manifold after the shift off of the manifold. Chapter 4. Predicting Voltage Instability of a Power System 43 t- o u.u E E" t- 1 Figure 4.4: View of reach tube at t = 0 (top) and t = 1 (bottom). In this image, the red surface (darker surface) is the manifold, the green surface (lighter curve) is the interface for the reachable tube and the blue surface represents the reachable tube. Chapter 4. Predicting Voltage Instability of a Power System 44 t- 2 t- 3 0^ 0.6, 0.4 v 0.2 v 0^  0.2.5 Figure 4.5: View of reach tube at t = 2 (top) and t = 3 (bottom). In this image, the red surface (darker surface) is the manifold, the green surface (lighter curve) is the interface for the reachable tube and the blue surface represents the reachable tube. Chapter 4. Predicting Voltage Instability of a Power System 45 1*4 Figure 4.6: View of reach tube at t = 4 (top) and t = 5 (bottom). In this image, the red surface (darker surface) is the manifold, the green surface (lighter curve) is the interface for the reachable tube and the blue surface represents the reachable tube. Chapter 4. Predicting Voltage Instability of a Power System 46 Chapter 4. Predicting Voltage Instability of a Power System 47 t- 5 Figure 4.8: The plot compares the interface for the reachable tube at the final timestep with simulations calculated by Matlab's D A E solver using the dynam-ics described in (4.1). The red surface represents the constraint manifold, the green surface represents the interface for the reachable tube and the blue lines represent the trajectories calculated by the D A E solver. Chapter 5 Conclusion 48 This thesis has made two contributions to the intellectual community. Firstly, it has added the fast marching method, signed distance and extension velocity functionality to Ian Mitchell's Level Set Toolbox. These functions use Dijkstra's algorithm to solve the non-linear Eikonal equation. This algorithm traverses each node in the grid once and it relies on a min-heap data structure to traverse these nodes in an optimal 0(N log N) time. The methods were implemented in Matlab in order to work with the Level Set Toolbox, but the min-heap struc-ture's operation time was very slow in this language so the functions had to be reimplemented in C and connected to Matlab via a M E X interface. This greatly improved the running time for all three methods. The implementation of the methods was proved accurate by running some well known example problems. Secondly, a method was developed, proven and tested for calculating reach-able sets on manifolds or equivalently DAEs. The algebraic constraint of the D A E represents the manifold and the differential constraint represents the dy-namics of the reachable set. A common practice for solving reachable sets in a full dimensional space is to represent the reachable set surface as an implicit surface function and to use level set methods to evolve the surface. The research in this thesis extends the level set computational technique by reinitializing the implicit surface function at each timestep with the closest point method.. This method naturally extends surface gradients off of the constraint manifold onto the entire computational domain, and it does so by assigning each point in the implicit surface function the same value as the value of the closest point on the constraint manifold. This procedure ensures that the gradients of the implicit surface function are constant in the normal direction of the constraint manifold. The closest point method was proven effective by solving two toy D A E problems and comparing the results to those calculated when the closest point method was not used. This method for solving reach sets on manifolds is used to solve a real life power generator example where the voltage dynamics are modeled with a D A E . This D A E has an additional problem that the model in inaccurate for a subset of the domain. In this subset the actual dynamics are index-2 or higher so the index-1 D A E which is being used to model the dynamics becomes singular at these points. Two fixes are done to the model to deal with the singular surface. The first fix damps out the dynamics inside the target set, since this set of points have already been marked as unsafe. The second fix uses the closest point method to extend the well modeled dynamics on the manifold throughout the entire domain; thus suggesting that the closest point method can be used Chapter 5. Conclusion 49 to extend both surface functions and velocity functions off of a manifold onto the entire domain in a manner which yields motion on the manifold consistent with the original dynamics. 50 Bibliography D. Adalsteinsson and J .A. Sethian. The fast construction of extension velocities in level set methods. Journal of Computational Physics, 148:2-22, 1999. Uri M . Ascher and Linda R. Petzold. Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations. Society for Industrial and Applied Mathematics, 1998. E . W. Dijkstra. A note on two problems in connection with graphs. Nu-merische Mathematik, 1:269-271, 1959. Ron Kimmel and James A . Sethian. Optimal algorithm for shape from shading and path planning. Journal of Mathematical Imaging and Vision, 14(3):224-241, 1992. Ian M . Mitchell. Comparing forward and backward reachability as tools for safety analysis. In Alberto Bemporad, Antonio Bicchi, and Giorgio But-tazzo, editors, Hybrid Systems: Computation and Control, number 4416, pages 428-443. Springer Verlag, 2007. Ian M . Mitchell. A toolbox of level set methods version 1.1. 2007. Stanley Osher and Ronald Fedkiw. Level Set Methods and Dynamic Implicit Surfaces. Springer-Verlag, 2003. J . Qian, Y . - T . Zhang, and H.-K. Zhao. Fast sweeping methods for eikonal equations on triangular meshes. SIAM Journal on Numerical Analysis, 45:83-107,2007. Elisabeth Rouy and Agnes Tourin. A viscosity solutions approach to shape-from-shading. SIAM Journal on Numerical Analysis, 29(3):867-884, June 1992. Steve J . Ruuth and Barry Merriman. A simple embedding method for solving partial differential equations on surfaces. 2006. http://www.math. s fu .ca /"sruuth/ . J .A. Sethian. A fast marching level set method for monotonically advancing fronts. Proceedings of the National Academy of Sciences, USA, 93(4): 1591-1595, 1996. Bibliography 51 [12] J .A. Sethian. Fast marching methods. SIAM Review, 41(2):199-235, July 1999. [13] Yoshihiko Susuki and Takashi Hikihara. Predicting voltage instability of power system via hybrid system reachability analysis. In Proceedings of the 2007 American Control Conference, pages 4166-4171, New York City, USA, July 2007. [14] Claire J . Tomlin, Ian Mitchell, Alexandre M . Bayen, and Meeko Oishi. Computational techniques for the verification of hybrid systems. Proceed-ings of the IEEE, 91:986-1001, July 2003. [15] Yen Hsi Richard Tsai. Rapid and accurate computation of the distance function using grids. Journal of Computational Physics, 178:175-195, 2002. [16] John N . Tsitsiklis. Efficient algorithms for globally optimal trajectories. IEEE Transactions on Automatic Control, 40(9):1528-1538, September 1995. • [17] V . Venkatasubramanian, H. Schattler, and J . Zaborszky. Voltage dynam-ics: Study of a generator with voltage control transmission, and matched mw load. IEEE Transactions on Automatic Control, 37(11):1717-1733, November 1992. 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            data-media="{[{embed.selectedMedia}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0052081/manifest

Comment

Related Items