1 static char help[] = "Trajectory sensitivity of a hybrid system with state-dependent switchings.\n";
2
3 /*
4 The dynamics is described by the ODE
5 u_t = A_i u
6
7 where A_1 = [ 1 -100
8 10 1 ],
9 A_2 = [ 1 10
10 -100 1 ].
11 The index i changes from 1 to 2 when u[1]=2.75u[0] and from 2 to 1 when u[1]=0.36u[0].
12 Initially u=[0 1]^T and i=1.
13
14 References:
15 + * - H. Zhang, S. Abhyankar, E. Constantinescu, M. Mihai, Discrete Adjoint Sensitivity Analysis of Hybrid Dynamical Systems With Switching,
16 IEEE Transactions on Circuits and Systems I: Regular Papers, 64(5), May 2017
17 - * - I. A. Hiskens, M.A. Pai, Trajectory Sensitivity Analysis of Hybrid Systems, IEEE Transactions on Circuits and Systems, Vol 47, No 2, February 2000
18 */
19
20 #include <petscts.h>
21
22 typedef struct {
23 PetscScalar lambda1;
24 PetscScalar lambda2;
25 PetscInt mode; /* mode flag*/
26 PetscReal print_time;
27 } AppCtx;
28
MyMonitor(TS ts,PetscInt stepnum,PetscReal time,Vec U,PetscCtx ctx)29 PetscErrorCode MyMonitor(TS ts, PetscInt stepnum, PetscReal time, Vec U, PetscCtx ctx)
30 {
31 AppCtx *actx = (AppCtx *)ctx;
32 Mat sp;
33 PetscScalar *u;
34 PetscInt nump;
35 FILE *f;
36
37 PetscFunctionBegin;
38 if (time >= actx->print_time) {
39 actx->print_time += 1. / 256.;
40 PetscCall(TSForwardGetSensitivities(ts, &nump, &sp));
41 PetscCall(MatDenseGetColumn(sp, 2, &u));
42 f = fopen("fwd_sp.out", "a");
43 PetscCall(PetscFPrintf(PETSC_COMM_WORLD, f, "%20.15lf %20.15lf %20.15lf\n", (double)time, (double)PetscRealPart(u[0]), (double)PetscRealPart(u[1])));
44 PetscCall(MatDenseRestoreColumn(sp, &u));
45 fclose(f);
46 }
47 PetscFunctionReturn(PETSC_SUCCESS);
48 }
49
EventFunction(TS ts,PetscReal t,Vec U,PetscReal * fvalue,PetscCtx ctx)50 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal *fvalue, PetscCtx ctx)
51 {
52 AppCtx *actx = (AppCtx *)ctx;
53 const PetscScalar *u;
54
55 PetscFunctionBegin;
56 PetscCall(VecGetArrayRead(U, &u));
57 if (actx->mode == 1) {
58 fvalue[0] = PetscRealPart(u[1] - actx->lambda1 * u[0]);
59 } else if (actx->mode == 2) {
60 fvalue[0] = PetscRealPart(u[1] - actx->lambda2 * u[0]);
61 }
62 PetscCall(VecRestoreArrayRead(U, &u));
63 PetscFunctionReturn(PETSC_SUCCESS);
64 }
65
ShiftGradients(TS ts,Vec U,AppCtx * actx)66 PetscErrorCode ShiftGradients(TS ts, Vec U, AppCtx *actx)
67 {
68 Mat sp;
69 PetscScalar *x;
70 PetscScalar *u;
71 PetscScalar tmp[2], A1[2][2], A2[2], denorm;
72 PetscInt nump;
73
74 PetscFunctionBegin;
75 PetscCall(TSForwardGetSensitivities(ts, &nump, &sp));
76 PetscCall(VecGetArray(U, &u));
77
78 if (actx->mode == 1) {
79 denorm = -actx->lambda1 * (u[0] - 100. * u[1]) + 1. * (10. * u[0] + u[1]);
80 A1[0][0] = 110. * u[1] * (-actx->lambda1) / denorm + 1.;
81 A1[1][0] = -110. * u[0] * (-actx->lambda1) / denorm;
82 A1[0][1] = 110. * u[1] * 1. / denorm;
83 A1[1][1] = -110. * u[0] * 1. / denorm + 1.;
84
85 A2[0] = 110. * u[1] * (-u[0]) / denorm;
86 A2[1] = -110. * u[0] * (-u[0]) / denorm;
87 } else {
88 denorm = -actx->lambda2 * (u[0] + 10. * u[1]) + 1. * (-100. * u[0] + u[1]);
89 A1[0][0] = 110. * u[1] * (actx->lambda2) / denorm + 1;
90 A1[1][0] = -110. * u[0] * (actx->lambda2) / denorm;
91 A1[0][1] = -110. * u[1] * 1. / denorm;
92 A1[1][1] = 110. * u[0] * 1. / denorm + 1.;
93
94 A2[0] = 0;
95 A2[1] = 0;
96 }
97
98 PetscCall(VecRestoreArray(U, &u));
99
100 PetscCall(MatDenseGetColumn(sp, 0, &x));
101 tmp[0] = A1[0][0] * x[0] + A1[0][1] * x[1];
102 tmp[1] = A1[1][0] * x[0] + A1[1][1] * x[1];
103 x[0] = tmp[0];
104 x[1] = tmp[1];
105 PetscCall(MatDenseRestoreColumn(sp, &x));
106
107 PetscCall(MatDenseGetColumn(sp, 1, &x));
108 tmp[0] = A1[0][0] * x[0] + A1[0][1] * x[1];
109 tmp[1] = A1[1][0] * x[0] + A1[1][1] * x[1];
110 x[0] = tmp[0];
111 x[1] = tmp[1];
112 PetscCall(MatDenseRestoreColumn(sp, &x));
113
114 PetscCall(MatDenseGetColumn(sp, 2, &x));
115 tmp[0] = A1[0][0] * x[0] + A1[0][1] * x[1];
116 tmp[1] = A1[1][0] * x[0] + A1[1][1] * x[1];
117 x[0] = tmp[0] + A2[0];
118 x[1] = tmp[1] + A2[1];
119 PetscCall(MatDenseRestoreColumn(sp, &x));
120 PetscFunctionReturn(PETSC_SUCCESS);
121 }
122
PostEventFunction(TS ts,PetscInt nevents,PetscInt event_list[],PetscReal t,Vec U,PetscBool forwardsolve,PetscCtx ctx)123 PetscErrorCode PostEventFunction(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec U, PetscBool forwardsolve, PetscCtx ctx)
124 {
125 AppCtx *actx = (AppCtx *)ctx;
126
127 PetscFunctionBegin;
128 /* PetscCall(VecView(U,PETSC_VIEWER_STDOUT_WORLD)); */
129 PetscCall(ShiftGradients(ts, U, actx));
130 if (actx->mode == 1) {
131 actx->mode = 2;
132 /* PetscCall(PetscPrintf(PETSC_COMM_SELF,"Change from mode 1 to 2 at t = %f \n",(double)t)); */
133 } else if (actx->mode == 2) {
134 actx->mode = 1;
135 /* PetscCall(PetscPrintf(PETSC_COMM_SELF,"Change from mode 2 to 1 at t = %f \n",(double)t)); */
136 }
137 PetscFunctionReturn(PETSC_SUCCESS);
138 }
139
140 /*
141 Defines the ODE passed to the ODE solver
142 */
IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,PetscCtx ctx)143 static PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, PetscCtx ctx)
144 {
145 AppCtx *actx = (AppCtx *)ctx;
146 PetscScalar *f;
147 const PetscScalar *u, *udot;
148
149 PetscFunctionBegin;
150 /* The next three lines allow us to access the entries of the vectors directly */
151 PetscCall(VecGetArrayRead(U, &u));
152 PetscCall(VecGetArrayRead(Udot, &udot));
153 PetscCall(VecGetArray(F, &f));
154
155 if (actx->mode == 1) {
156 f[0] = udot[0] - u[0] + 100 * u[1];
157 f[1] = udot[1] - 10 * u[0] - u[1];
158 } else if (actx->mode == 2) {
159 f[0] = udot[0] - u[0] - 10 * u[1];
160 f[1] = udot[1] + 100 * u[0] - u[1];
161 }
162
163 PetscCall(VecRestoreArrayRead(U, &u));
164 PetscCall(VecRestoreArrayRead(Udot, &udot));
165 PetscCall(VecRestoreArray(F, &f));
166 PetscFunctionReturn(PETSC_SUCCESS);
167 }
168
169 /*
170 Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
171 */
IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,PetscCtx ctx)172 static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, PetscCtx ctx)
173 {
174 AppCtx *actx = (AppCtx *)ctx;
175 PetscInt rowcol[] = {0, 1};
176 PetscScalar J[2][2];
177 const PetscScalar *u, *udot;
178
179 PetscFunctionBegin;
180 PetscCall(VecGetArrayRead(U, &u));
181 PetscCall(VecGetArrayRead(Udot, &udot));
182
183 if (actx->mode == 1) {
184 J[0][0] = a - 1;
185 J[0][1] = 100;
186 J[1][0] = -10;
187 J[1][1] = a - 1;
188 } else if (actx->mode == 2) {
189 J[0][0] = a - 1;
190 J[0][1] = -10;
191 J[1][0] = 100;
192 J[1][1] = a - 1;
193 }
194 PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));
195
196 PetscCall(VecRestoreArrayRead(U, &u));
197 PetscCall(VecRestoreArrayRead(Udot, &udot));
198
199 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
200 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
201 if (A != B) {
202 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
203 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
204 }
205 PetscFunctionReturn(PETSC_SUCCESS);
206 }
207
208 /* Matrix JacobianP is constant so that it only needs to be evaluated once */
RHSJacobianP(TS ts,PetscReal t,Vec X,Mat Ap,PetscCtx ctx)209 static PetscErrorCode RHSJacobianP(TS ts, PetscReal t, Vec X, Mat Ap, PetscCtx ctx)
210 {
211 PetscFunctionBeginUser;
212 PetscFunctionReturn(PETSC_SUCCESS);
213 }
214
main(int argc,char ** argv)215 int main(int argc, char **argv)
216 {
217 TS ts; /* ODE integrator */
218 Vec U; /* solution will be stored here */
219 Mat A; /* Jacobian matrix */
220 Mat Ap; /* Jacobian dfdp */
221 PetscMPIInt size;
222 PetscInt n = 2;
223 PetscScalar *u;
224 AppCtx app;
225 PetscInt direction[1];
226 PetscBool terminate[1];
227 Mat sp;
228 PetscReal tend;
229 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
230 Initialize program
231 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
232 PetscFunctionBeginUser;
233 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
234 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
235 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
236 app.mode = 1;
237 app.lambda1 = 2.75;
238 app.lambda2 = 0.36;
239 app.print_time = 1. / 256.;
240 tend = 0.125;
241 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex1fwd options", "");
242 {
243 PetscCall(PetscOptionsReal("-lambda1", "", "", app.lambda1, &app.lambda1, NULL));
244 PetscCall(PetscOptionsReal("-lambda2", "", "", app.lambda2, &app.lambda2, NULL));
245 PetscCall(PetscOptionsReal("-tend", "", "", tend, &tend, NULL));
246 }
247 PetscOptionsEnd();
248
249 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
250 Create necessary matrix and vectors
251 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
252 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
253 PetscCall(MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
254 PetscCall(MatSetType(A, MATDENSE));
255 PetscCall(MatSetFromOptions(A));
256 PetscCall(MatSetUp(A));
257
258 PetscCall(MatCreateVecs(A, &U, NULL));
259
260 PetscCall(MatCreate(PETSC_COMM_WORLD, &Ap));
261 PetscCall(MatSetSizes(Ap, n, 3, PETSC_DETERMINE, PETSC_DETERMINE));
262 PetscCall(MatSetType(Ap, MATDENSE));
263 PetscCall(MatSetFromOptions(Ap));
264 PetscCall(MatSetUp(Ap));
265 PetscCall(MatZeroEntries(Ap));
266
267 PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, n, 3, NULL, &sp));
268 PetscCall(MatZeroEntries(sp));
269 PetscCall(MatShift(sp, 1.0));
270
271 PetscCall(VecGetArray(U, &u));
272 u[0] = 0;
273 u[1] = 1;
274 PetscCall(VecRestoreArray(U, &u));
275
276 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
277 Create timestepping solver context
278 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
279 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
280 PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
281 PetscCall(TSSetType(ts, TSCN));
282 PetscCall(TSSetIFunction(ts, NULL, (TSIFunctionFn *)IFunction, &app));
283 PetscCall(TSSetIJacobian(ts, A, A, (TSIJacobianFn *)IJacobian, &app));
284
285 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
286 Set initial conditions
287 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
288 PetscCall(TSSetSolution(ts, U));
289
290 PetscCall(TSForwardSetSensitivities(ts, 3, sp));
291 /* Set RHS JacobianP */
292 PetscCall(TSSetRHSJacobianP(ts, Ap, RHSJacobianP, &app));
293
294 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
295 Set solver options
296 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
297 PetscCall(TSSetMaxTime(ts, tend));
298 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
299 PetscCall(TSSetTimeStep(ts, 1. / 256.));
300 PetscCall(TSMonitorSet(ts, MyMonitor, &app, NULL));
301 PetscCall(TSSetFromOptions(ts));
302
303 /* Set directions and terminate flags for the two events */
304 direction[0] = 0;
305 terminate[0] = PETSC_FALSE;
306 PetscCall(TSSetEventHandler(ts, 1, direction, terminate, EventFunction, PostEventFunction, (void *)&app));
307 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
308 Run timestepping solver
309 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
310 PetscCall(TSSolve(ts, U));
311
312 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
313 Free work space. All PETSc objects should be destroyed when they are no longer needed.
314 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
315 PetscCall(MatDestroy(&A));
316 PetscCall(VecDestroy(&U));
317 PetscCall(TSDestroy(&ts));
318
319 PetscCall(MatDestroy(&Ap));
320 PetscCall(MatDestroy(&sp));
321 PetscCall(PetscFinalize());
322 return 0;
323 }
324
325 /*TEST
326
327 build:
328 requires: !complex
329
330 test:
331 args: -ts_monitor
332
333 TEST*/
334