User's guide for the ROC-HJ solver: Finite Differences and Semi
Contents
1. e function compute_Hconst This function must define d constants which are bounds for ot a p i i 1 d for x Q and for possible gradient values p and time t This is then used to set the time step At in order to satisfy a CFL condition These constants may also be used in the function Hnum The user may also define directly the constants c Hconst i and initialize the previous function accordingly Then in all cases we have also to set the scheme discretisation parameters e TYPE SCHEME LF ENO2 type of spacial discretization LF Lax Friedrich scheme first order scheme ENO2 ENO scheme of second order to approximate the derivatives Vu e TYPE RK RK1 RK2 RK3 Time discretization by a Runge Kutta method of order 1 2 or 3 RK1 RK method of order 1 RK2 RK method of order 2 RK3 RK method of order 3 Semi Lagrangian Method This method is used when METHOD MSL It assumes that H a2 Vu max f z a Vu L x a Then it needs to define e function dynamics dynamics f z a e function distributed cost cost function L x a e TYPE_SCHEME STA EVO STA stationnary case for solving 5 In this case the scheme is based on an iterative procedure EVO dynamic case for solving 1 ORDER this value is set to 1 for solving first order equations of type 1 see other section below for second order HJB equations e TYPE STA LOOP NORMAL SPECIAL Mesh loop2. use exe nt 4 3 User s parameters The file data data_simulation h is the main user s input file and contains the parame ters that define the equation to be solved It can include an other data xxx h file see basic examples data_basicmodel h data_FD_2d_ex1_basic h It contains also parame ters relative to the output during the execution saving data and data post treatment see paragraph Output parameters e NAME name of the problem e DIM dimension d of the problem Choice of the solver e METHOD MDF MSL MDF Finite Difference Method MSL Semi Lagrangian Method Stopping criteria e T terminal time e MAX_ITERATION to stop the program when this maximum number of iterations is reached e EPSILON for stopping criteria If set to 0 do nothing If set to some positive value then the program will stop at iteration n 1 as soon as utt p00 max v yP lt e 2 the LY error between two successive steps is smaller than Finite Difference Method This method is used when METHOD MDF and it utilizes e COMMANDS 0 1 2 This decides how we define the Hamiltonian function H x p using the commands or not 1 the commands are used in the definition of the numerical hamiltonian as in 2 The program will use a numerical hamiltonian function corresponding to a finite difference approximation of H x Vu max f z a Vu x a Examples are given in data_basicmodel
3. will be present in the file OUTPUT VEX dat contains the final value as programmed in Vex if COMPUTE_VEX 1 OUTPUT coupe dat and coupeex dat can be obtained in the same way for a pro blem in dimension d gt 3 and when we make a cut in some particular subspace with some given coordinates and that we want to visualize the results See para graph Output parameters for using COUPE_DIMS and COUPE_VAL This is typically used to make a 2d cut see for instance example data data_FD_3d h 6 if V tn x is the value function at time tn attemp to reconstruct a trajectory such that infa V tn41 Y 2 tn41 is minimal this in the case OPTIM MAXIMUM 10 Furthermore the following files are related to trajectory reconstruction e QUTPUT traj n dat in case TRAJPT gt 1 corresponds to the trajectory number n Each lines has the form XI X2 sas R al vos ap tv where xi are the coordinates of the point at time t and aj the corresponding control parameters e OUTPUT successTrajectories dat Each line n of this file has the form xi x2 xd t success and indicates weither the trajectory number n was successfull or not success 1 or 0 and shows the corresponding terminal coordinates xi and time t of the last constructed point 5 Graphic outputs Matlab the file OUTPUT output_view m can be executed by typing output view in the Matlab s user interface we advise to first set the current Matlab s directory
4. g t 2 0 6 where g t x is defined by the user In particular it forces to have u t x gt g t x in the case of 6 Also the following obstacle equations can be considered max amp H x Vu u ata 0 7 max min F H x Vu u u ts a att 0 8 where g t x is used defined We will have u t x lt g t x in the case of 7 or u t x g t x g t x in the case of 8 For this the following parameters are used e OBSTACLE 0 1 0 no obstacle g taken into account 1 Equation 6 is treated with the obstacle g e function g_obstacle e OBSTACLE_TILDE 0 1 0 no obstacle g taken into account 1 Equation 7 is treated with the obstacle g e function g_obstacle_tilde 4 Or OUTPUT VF_PROCxx dat values if parallel MPI is used In the case we set OBSTACLE 1 and OBSTACLE_TILDE 1 then 8 is treated and both functions g and g must be defined The obstacle functions should satisfy g x lt g a in order to avoid undesidered results Output parameters output files are generaly contained in the directory OUTPUT COMPUTE_MAIN_LOOP 1 the main computational loop iterative scheme is called 0 the main loop is not called only data initializations and savings are performed COMPUTE_VEXE 0 1 will save a VEX dat file 1 if the exact solution is given by the user using function Vex 0 no saving In particular it will compute errors see below relative
5. f o x a r t z aju b t x a Vu sTr ott x a jo t z at Deu ae te 0 7 z R u 0 x uo z TER 4b where A is some non empty compact subset of R m gt 1 of the form 4 ai il b t x a is a vector of RY r t x a and L t x a are real valued and a t x a is ad xp real matrix for some p gt 1 This problem is linked to the computation of the value function of stochastic optimal control problems 1 2 Steady equations The problem is to find u u x solution of A x u a H z u x Vu 2 0 TEN 5 with H given directly or in the form of 2 or as in 4a for second order equations the function A x should be strictly positive Q an hyperrectangle of R of the form TZ a bil Equation 5 is complemented by boundary conditions as follows u x Joora X x E ON It is possible also to solve 5 in a subdomain C In that case the set C should be defined such that C gaomain x lt 0 equivalently Q C x Q gJdomain gt 0 and the boundary conditions at the border OC should be of Dirichlet type fixed by uo x u x uox xe OC In general this equation is solved by using an iterative procedure value iteration where the part x u is treated explicitly for a FD method or implicitly for a SL method 3 2 MaXae4 Minpew can be also minac a Maxyep 3 For steady equations the scheme is based on the following iteration on n gt
6. h Goss or data_FD_2d_ex1_basic h eikonal equation 2 the commands are used in the definition of the numerical hamiltonian as in 3 The program will use a numerical hamiltonian function corresponding to a finite difference approximation of for instance H z Vu max min f s a B Vu l x a B 0 an Hnum function must be directly defined by the user it has to be a numerical hamiltonian corresponding to H z p Examples are given in data_advancedmodel data_FD_2d_exi_advanced h If COMMANDS 1 we need also to define e function dynamics f z a e function distributed cost L x a e aset of commands corresponding to the discretization of a cube heed lai Bi defi ned by cDIM dimension ne of the set of controls NCD cDIM number of commands per direction UMIN cDIM UMAX cDIM min and max values q and 3 of the com mands in each direction If COMMANDS 2 we need to define e function dynamics2 f zx a b t e function distributed cost2 amp x a b t e two sets of commands A and B corresponding to the discretization of A rarer o4 BA and B thas aP BPI The parameters cDIM cDIM2 will correspond to n4 and np the dimension of each set of controls NCD cDIM is the number of commands per direction and UMIN cDIM UMAX cDIM are the min and max values of a and Bi of the commands in each direction If COMMANDS 0 we need to define e function Hnum a first order numerical Hamiltonian
7. to OUTPUT aome parameters are given below For the first figure e level set a level set value e PLOT_CONTOURE 0 1 if set to 1 will plot contours of the output with level set value e PLOT_REACHABLEE 0 1 if set to 1 will plot the 2d set of points where output value is lt level_set This plot is based on the output value which is in the file VF dat for DIM 2 and a priori in the file coupe dat for DIM gt 2 For the second figure e PLOT 3d 10 1 set this parameter to 1 in order to obtain a 3d graph of the result e PLOT TMINE 0 1 default is 0 Set this parameter to 1 in order to obtain a 3d graph of the minimal time to reach a target Will then use output tmin dat and will draw several contour plots ranging from 0 to T plotting tmin dat is only ok for DIM 2 Also for both figures e AXIS_EQUALE 0 1 if set to 1 will make axis equal e TRAJECTORYE 0 1 if set to 1 will line plot the list of first two components x1 2 of file trajectory dat 11 Paraview For the moment this option is only working for some 2d 3d cases Launch paraview then load the files in VTK load tab and so on 6 Examples Rotation example We give two ways to program an advection like example in files data_basicmodel h and data_advancedmodel h We consider the equation Ou max 0 f x Vu 0 with the parameters Q 2 2 T 0 5 f a1 2 27 a2 x1 Hence the Hamilto ni
8. 0 for a fixed At gt 0 urt gx u x At Alxju x h x u x Du x 0 2 C 2 Compilation and execution under GCC and CMake In order to make the program works the c compiler GCC and the build system CMake are proposed to help the compilation process The source files are typically in the directory src User s data files are typically in the directory data For deeper informations we redirect the user to the following links http gcec gnu org install http www cmake org cmake resources software html Building the Makefile cmake do not forget the dot Compilation make source cpp files are generally in the directory src Options make clean cleaning some o files make cleanall Other clean to first erase all unecessary files can be done before the cmake command Sequential code Execution basic exe Execution with options exe nn NN nc NC Options e option nn NN number of mesh points per dimension e option nc NC number of controls per command s dimension where h x u x Du x corresponds to some numerical approximation of H x u x Vu x Iterations are stopped when u u z is smaller than a given threshold Therefore we obtain the following recursion u t x u x At eura h x u x Du 2 0 EC For the SL method the iterative scheme is based on an implicit treatment of the A
9. User s guide for the ROC HJ solver Finite Differences and Semi Lagrangian methods November 2014 Version 2 07 Current developers O Bokanowski H Zidani J Zhao A Desilles 1 Problem solved This is a c MPI OpenMP library for solving d dimensional Hamilton Jacobi Bellman equations by finite difference methods or semi lagrangian methods First order and second order HJ equations time dependent or steady equations can be solved More precisely the equations considered are of the following type 1 1 Evolutionary case The problem is to find u u t x solution of ou X x u x H t u x Vu 0 TEN te 0 7 u t oralt x EON te 0 7 1 u 0 x u x TERN where Q is a domain of R of the form Tesla bi and gbora and uo are given functions other boundary conditions are possible The function H t 2 p can be defined either directly or as follows t H t x Vu max f t 2 a Vu L t 2 i 2 acA Reachability Optimal Control and Hamilton Jacobi equations tIncluding a SL solver for second order HJ equations see Section 7 1 maxge can be also minge4 where A is a set of controls of the form 4 ai 8i or A t x Vu max min f t x a b Vu L t x a D 3 acA beB where A and B are control sets of the form J 4 a 64 and i2 la BF Some second order equations can also be treated of the following type Ou 0 En A x u 4a max
10. a is a d x p real matrix for some p gt 1 This problem is linked to the computation of the value function of stochastic optimal control problems It is also possible to consider the corresponding steady equation that is equation 12a i Ou alone with no term ot The c proposed solver is based on an SL method 7 User inputs ORDER 2 PARAMP denoted p below This is in general set to p d in the scheme function Sigma X x k a t is the kth column vector of the matrix 4 t x a function Drift B x k a t corresponding to a set a set of p vectors B t x a k 1 p of R such that a B t 2 a B t x a k 1 p For instance the user can set B t 2 a B t z a and B 0 Vk gt 2 or BE the k th component of B the other components beeing set to zero if p d e function funcR r t x a e function distributed cost x t x distributive cost L t x a e function funcL x A x assuming A x gt r t x a in the case of a steady equation that is equation 12a alone with no term ou Scheme definition we consider the following SL scheme implemented on the points x of the grid G The initialization is done by v 7 uo x EG 7 For advanced users it is programmed in function secondorder_itSL_evo that might be available in the source file HJIB_SL cpp or in include secondorder_itSL_evo h 13 For n 0 Nr 1 time iterations or untill some st
11. an H z p t max 0 f p In the first way data_basicmodel h we define dynamics f a a af x with two control values a 0 1 In the second way data_advancedmodel h we define the Lax Friedrich numerical hamiltonian Hnum associated to H see 11 Eikonal equation seethe files data_FD_2d_ex1_basic h and data_FD_2d_ex1_advanced h The problem solved is u c x t Vul 0 a 2 2 t 0 7 9 u 0 x uox x 2 21 10 with d 2 T 1 here c x t 1 and with some radially symetric initial data uo x In the file data_FD_2d_ex1_basic h a scheme is programmed using a control discretisation of Vul as follows 2kr Vul aia aS Conin sin Vah 6 fe where NCD is the number of controls In the file data FD 2d ex1 advanced h a scheme is programmed using a Lax Friedriech numerical approximation Hnum associated to H _ z p_ pt Pi Pi Hrum t PI pT OP H z Sa S 5 11 1 1 The exact solution is given for comparison 12 7 Second order HJB equations and SL scheme The problem solved is the following second order Hamilton Jacobi HJ equation Ou FE A x u 12a 1 max 5Tr X t x a t x 4 D u B t x a Vu r t x a u L t x a ae te 0 7 x RI 0 r mir TER 12b where A is some non empty compact subset of R m gt 1 B t x a is a vector of RY r t xz a amp t x a are real valued and X t x
12. is 2 NC total number of commands NCD cDIM tab indicating the number of discretisation points for each control u i UMIN UMAX bounds for each control parameter Then the following more specific parameters are used TRAJPT 0 1 set to 1 for trajectory reconstruction 0 otherwise initialpoint DIM coordinate of initial point xo for trajectory reconstruction TRAJ_METHOD method of reconstruction 0 if based on minimal time utilizes tmin dat 1 if based on the value function utilizes VFxx dat files time_TRAJ_START starting time to used in the case TRAJ_METHOD 1 Can be any value in 0 T The program will construct an approximated trajectory such that y t f t y t a t t to T and starting with y to zo as initial point Note that several trajectory reconstructions can be obtained by setting TRAJPT gt 1 equal to the desired number of trajectories and then initialpoint TRAJPT DIM should contain list of corresponding points If TRAPJT 0 then one should define accordingly initialpoint DIM Output files Output files are generaly contained in the directory OUTPUT OUTPUT VF dat contains the final u at the end of the computation The file is struc tured as follows on each line case FORMAT_FULLDATA 1 il i2 id val where val corresponds to the value u T x at mesh point x j Liz In the case FORMAT_FULLDATA 0 on each line in the same order only the value val
13. ly to this value COMPUTE_TMINE 0 1 will compute file tmin dat which contains the first time tn 0 T or some P1 interpolant such that v tn x lt 0 Set the value to INF 1 e5 otherwise PRECOMPUTE_COORD 1 precomputes the coordinates faster computations but more memory deman ding 0 no precomputations SAVE_VFALL 0 1 to save the value of v each SAVE_VFALL_STEP iterations into files VFxx dat FORMAT_FULLDATA 0 1 if set to 0 to save only the value at grid mesh points for files such as VF dat CHECK ERROR 0 1 to compute error every CHECK ERROR STEP steps LY and Lt error computations Errors are relative to the Vex function Error are computed only in the region of points x such that Vex t x lt C_THRESHOLD COUPE_DIMS DIM COUPE_VALS DIM list of integers n 0 or 1 to define the two variables used for the cut list of values c or 0 to define the position of the cut The values c are used only when n 0 The output is in the file output dat Example COUPE_DIMS DIM 1 0 1 and COUPE_VALS DIM 0 0 5 0 define a cut in the plane x2 0 5 Parameters for trajectory reconstruction The dynamics used is the one defined in dynamics 5 For the moment C_THRESHOLD is defined in src main cpp The following parameters are for control discretization they are similar to the one used in the case COMMANDS 1 or COMMANDS 2 cDIM dimension of the control default
14. opping criteria is satisfied in the case of steady equations for all grid points x G we consider unti eT txa At 1 2 acd 1 A a At XO ol ybe At At ee x a 13 2p o R Pp where the characteristics y h can be defined at iteration t tn for some h gt 0 for instance by yE h x BF t x a h eD t 2 a Vh 1 1 k 1 p 14 Also v y denotes some interpolation of v at point y typically Q1 8 8 Source architecture In this section we present the architecture used in the implementation of the pro ject The HJB class represents the main class and the DF and SL classes are sub classes of HJB Commands and Mesh are seperated classes We remind that the user needs to code the problem parameters into a C header file and to include its name in the file data data_simulation h For example if the user data file is called data_myexample h then the line include data_myexample h and only this one must be included in the file data data_simulation h 8 When using the definition 14 the scheme is of expected order O At O where e is of the order of the interpolation error v v For smooth data this interpolation error is or order Ax for Q interpolation where Az is the spatial mesh size 14 data h A include data h User file 15
15. or the time step For the semi lagrangian approach and for the stationnary equation the parameter h in the iterative procedure is also set to DT Other parameters for the definition of the problem The parameters are e XMIN DIM XMAX DIM the boundary a b of the domain in each direction e PERIODIC DIM each componant set the periodicity for each direction 1 direction is periodic O direction is not periodic e BOUNDARY 0 1 2 If set to 0 no boundary treatment If set to 1 Utilizes Dirichlet boundary conditions for each direction d such that PERIODIC d 0 The boundary value is given in the user defined function g bord If set to 2 Utilizes a special Ur zx 0 boundary condition for each direction d such that PERIODIC d 0 note that this special approximation can be unstable e function vo the initial data e function Vex if known the user can put the exact solution of the problem here e integer EXTERNAL vo 0 1 if O initialize with function vo otherwise initialize with last VF dat values located in OUTPUT VF dat 4 e COMPUTE_IN_SUBDOMAIN 0 1 determines if subdomain computations should be done so that the evaluation of u t x is done only at grid mesh point x such that Jdomain X lt 0 e function g_domain x function used to define the subdomain Obstacle terms Instead of solving u H t x u Vu 0 the solver can also treat HJ obstacle equations such as min 5 H t x u Vu u
16. order during mainloop for the stationary case NORMAL normal ordering loop SPECIAL special ordering loop makes 2 loops at each iteation modifiying the current values of the data v during each loop e TYPE RK RK1_EULER RK2_HEUN RK2_PM RK1_EULER RK1 Euler scheme RK2_HEUN RK2 Heun scheme RK2 PM RK2 Mid point scheme e INTERPOLATION BILINEAR PRECOMPBL DIRPERDIR BILINEAR default value bilinear interpolation Q1 interpolation PRECOMPBL precompute the interpolation coefficients to use with tiny mesh sizes DIRPERDIR another interpolation method e COMPLEMENTS P_INTERMEDIATE number of intermediate timesteps to go from tn to tn At for computing a caracteristic for a given RK method ORDER 1 2 decide treatement for first order HJ equation or for second order HJ equation Discretization parameters e ND DIM tab containing the size mesh in each direction cartesian mesh e MESH 0 1 default is 1 It utilizes ND i 1 points in direction i If MESH 0 the mesh points are at the center of the mesh cells and there are ND i points in direction i e DT the time step used for the solver for the evolutive equation For the finite difference approach if DT 0 then time step DT is computed so that the CFL condition be satisfied that is such that DT Hconst 0 dx 0 Hconst 1 dx 1 lt CFL where the CFL number belongs to 0 1 Otherwise if DT gt 0 then the value DT is used f
17. x u x term and becomes in the case of 2 urt 7 min m s wrie Atf t x a ut 2 a 0 ref Here u denotes the Q interpolation of u on the spatial grid mesh OpenMP Some options run with OpenMP Execution exe nt nbth where nbth is the number of threads typically try nbth 2 or nbth 4 Parallel code Parallel MPI version only Execution mpirun n MPI_PROCS exe nn NN nc NC nt OMP_THREADS nd MPI_GRID_DIM where e MPI_PROCS number of MPI processors default 1 e NN number of mesh points per dimension e NC number of controls per command s dimension e OMP_THREADS number of OpenMP threads default 1 e MPI_GRID_DIM dimension of the MPI mesh grid decomposition 2 or 3 default 2 Example 1 mpirun n 2 exe nn 200 will execute the program with 2 MPI processors and 200 mesh points per dimension Example 2 mpirun n 64 exe nn 500 nc 10 nt 4 nd 2 will execute the pro gram with MPI_PROCS 64 MPI processors NN 500 points per dimension NC 10 controls per command s dimension OMP_THREADS 4 OpenMP threads and MPI_GRID_DIM 2 the dimension of the MPI mesh grid decomposition Example 3 mpirun n 1 exe nn 500 nc 10 nt 4 The OpenMP parallelization is working for Semi Lagrangian and Finite Difference me thods When using only this method one can indicate only one process for the mpirun command for instance mpirun n 1 exe nt 4 for using 4 threads or equivalently
Download Pdf Manuals
Related Search