From abd7342d238a852eda9875383ec198da2b635ce2 Mon Sep 17 00:00:00 2001 From: "Documenter.jl" Date: Wed, 14 Aug 2024 13:30:06 +0000 Subject: [PATCH] build based on 1655a9d --- dev/.documenter-siteinfo.json | 2 +- ...e-28f28577.svg => 2D-example-58013bc8.svg} | 192 +- ...e-506e034f.svg => 2D-example-99a24f1e.svg} | 192 +- ...e-37751a20.svg => 2D-example-e4390f57.svg} | 192 +- ...e-0cd43820.svg => 2D-example-f2116a06.svg} | 58 +- dev/2D-example.html | 140 +- ...978.svg => 2D-preconditioner-6038766b.svg} | 4170 ++++++++--------- ...e80.svg => 2D-preconditioner-a1c3d910.svg} | 328 +- ...d1f.svg => 2D-preconditioner-d0a533dd.svg} | 320 +- ...bd5.svg => 2D-preconditioner-d454d92a.svg} | 254 +- dev/2D-preconditioner.html | 82 +- dev/index.html | 4 +- dev/search_index.js | 2 +- 13 files changed, 2968 insertions(+), 2968 deletions(-) rename dev/{2D-example-28f28577.svg => 2D-example-58013bc8.svg} (80%) rename dev/{2D-example-506e034f.svg => 2D-example-99a24f1e.svg} (86%) rename dev/{2D-example-37751a20.svg => 2D-example-e4390f57.svg} (80%) rename dev/{2D-example-0cd43820.svg => 2D-example-f2116a06.svg} (87%) rename dev/{2D-preconditioner-246e2978.svg => 2D-preconditioner-6038766b.svg} (64%) rename dev/{2D-preconditioner-2ee2be80.svg => 2D-preconditioner-a1c3d910.svg} (86%) rename dev/{2D-preconditioner-e66ead1f.svg => 2D-preconditioner-d0a533dd.svg} (85%) rename dev/{2D-preconditioner-b8e7ebd5.svg => 2D-preconditioner-d454d92a.svg} (81%) diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json index bb5dbb8..2575f8c 100644 --- a/dev/.documenter-siteinfo.json +++ b/dev/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.10.4","generation_timestamp":"2024-08-14T07:40:28","documenter_version":"1.5.0"}} \ No newline at end of file +{"documenter":{"julia_version":"1.10.4","generation_timestamp":"2024-08-14T13:30:02","documenter_version":"1.5.0"}} \ No newline at end of file diff --git a/dev/2D-example-28f28577.svg b/dev/2D-example-58013bc8.svg similarity index 80% rename from dev/2D-example-28f28577.svg rename to dev/2D-example-58013bc8.svg index c24236c..5079091 100644 --- a/dev/2D-example-28f28577.svg +++ b/dev/2D-example-58013bc8.svg @@ -1,118 +1,118 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/2D-example-506e034f.svg b/dev/2D-example-99a24f1e.svg similarity index 86% rename from dev/2D-example-506e034f.svg rename to dev/2D-example-99a24f1e.svg index 91fc0b2..897b0b0 100644 --- a/dev/2D-example-506e034f.svg +++ b/dev/2D-example-99a24f1e.svg @@ -1,118 +1,118 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/2D-example-37751a20.svg b/dev/2D-example-e4390f57.svg similarity index 80% rename from dev/2D-example-37751a20.svg rename to dev/2D-example-e4390f57.svg index ff6c37d..304eb5e 100644 --- a/dev/2D-example-37751a20.svg +++ b/dev/2D-example-e4390f57.svg @@ -1,118 +1,118 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/2D-example-0cd43820.svg b/dev/2D-example-f2116a06.svg similarity index 87% rename from dev/2D-example-0cd43820.svg rename to dev/2D-example-f2116a06.svg index 4fb9277..e3fedf2 100644 --- a/dev/2D-example-0cd43820.svg +++ b/dev/2D-example-f2116a06.svg @@ -1,44 +1,44 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/2D-example.html b/dev/2D-example.html index 2c03071..63e34e4 100644 --- a/dev/2D-example.html +++ b/dev/2D-example.html @@ -26,45 +26,45 @@ # Plot plot(range(-7, 2, 500), S, xlim = [-7, 2]) plot!([-7,2], [0,0], color = :black) -plot!(xlabel = "p0", ylabel = "S", legend=false)Example block output

Finite difference method

The main goal now is to find the zero of $S$. To this purpose, we use the numerical solver $\texttt{hybrd1}$ given in the package $\texttt{MINPACK.jl}$. If we don't provide the Jacobian $J_S$ of $S$ to the solver, the finite difference method is used to approximate it. 

ξ = [-1.0]                                                  # initial guess
+plot!(xlabel = "p0", ylabel = "S", legend=false)
Example block output

Finite difference method

The main goal now is to find the zero of $S$. To this purpose, we use the numerical solver $\texttt{hybrd1}$ given in the package $\texttt{MINPACK.jl}$. If we don't provide the Jacobian $J_S$ of $S$ to the solver, the finite difference method is used to approximate it. 

ξ = [-1.0]                                                  # initial guess
 S!(s, ξ) = (s[:] .= S(ξ[1]); nothing)                       # intermediate function
 p0_sol = fsolve(S!, ξ, show_trace = true)                   # solve
 println(p0_sol)
Iter     f(x) inf-norm    Step 2-norm      Step time
 ------   --------------   --------------   --------------
-     1     3.000000e+00     0.000000e+00         0.167741
-     2     9.318967e-01     1.069263e+00         0.224512
-     3     3.908179e-08     2.171079e-01         0.001074
-     4     3.126006e-08     3.818465e-16         0.000962
-     5     2.602624e-08     7.541206e-17         0.000934
-     6     1.122729e-07     1.864781e-15         0.000952
-     7     1.251064e-08     1.228963e-15         0.000947
-     8     1.071865e-08     6.960129e-18         0.000969
-     9     3.372047e-08     2.490162e-16         0.000917
-    10     6.149611e-09     1.433785e-16         0.000945
-    11     2.836331e-08     2.624355e-17         0.000894
-    12     2.550366e-08     1.772446e-17         0.000881
-    13     2.794534e-08     1.292187e-17         0.001759
-    14     2.543587e-08     1.364932e-17         0.000882
-    15     6.042093e-09     4.286522e-19         0.000957
-    16     3.036848e-08     7.912620e-17         0.000898
-    17     2.599380e-08     5.504418e-17         0.000906
-    18     2.518026e-08     1.434518e-18         0.000887
-    19     6.005500e-09     5.041160e-20         0.000960
-    20     3.036868e-08     7.817582e-17         0.000889
-    21     2.601931e-08     5.449264e-17         0.000909
-    22     2.521370e-08     1.406704e-18         0.000919
-    23     5.969729e-09     4.888090e-20         0.018251
-    24     3.036827e-08     7.723669e-17         0.001152
-    25     2.604422e-08     5.394382e-17         0.000918
-    26     2.524639e-08     1.379606e-18         0.000882
-    27     5.934758e-09     4.741150e-20         0.000949
-    28     3.036860e-08     7.634290e-17         0.000919
-    29     2.606871e-08     5.342255e-17         0.000911
-    30     2.527841e-08     1.353711e-18         0.000896
-    31     5.900549e-09     4.601691e-20         0.000998
-    32     3.036849e-08     7.546245e-17         0.000948
-    33     2.609265e-08     5.290603e-17         0.000948
-    34     2.530975e-08     1.328510e-18         0.000889
+     1     3.000000e+00     0.000000e+00         0.191782
+     2     9.318967e-01     1.069263e+00         0.235390
+     3     3.908179e-08     2.171079e-01         0.001081
+     4     3.126006e-08     3.818465e-16         0.000963
+     5     2.602624e-08     7.541206e-17         0.000942
+     6     1.122729e-07     1.864781e-15         0.000981
+     7     1.251064e-08     1.228963e-15         0.000962
+     8     1.071865e-08     6.960129e-18         0.000954
+     9     3.372047e-08     2.490162e-16         0.000922
+    10     6.149611e-09     1.433785e-16         0.000987
+    11     2.836331e-08     2.624355e-17         0.000917
+    12     2.550366e-08     1.772446e-17         0.000969
+    13     2.794534e-08     1.292187e-17         0.001826
+    14     2.543587e-08     1.364932e-17         0.000933
+    15     6.042093e-09     4.286522e-19         0.000989
+    16     3.036848e-08     7.912620e-17         0.000944
+    17     2.599380e-08     5.504418e-17         0.000942
+    18     2.518026e-08     1.434518e-18         0.000919
+    19     6.005500e-09     5.041160e-20         0.000969
+    20     3.036868e-08     7.817582e-17         0.000913
+    21     2.601931e-08     5.449264e-17         0.000922
+    22     2.521370e-08     1.406704e-18         0.000895
+    23     5.969729e-09     4.888090e-20         0.000982
+    24     3.036827e-08     7.723669e-17         0.000931
+    25     2.604422e-08     5.394382e-17         0.000914
+    26     2.524639e-08     1.379606e-18         0.000905
+    27     5.934758e-09     4.741150e-20         0.000974
+    28     3.036860e-08     7.634290e-17         0.000893
+    29     2.606871e-08     5.342255e-17         0.000930
+    30     2.527841e-08     1.353711e-18         0.000905
+    31     5.900549e-09     4.601691e-20         0.000974
+    32     3.036849e-08     7.546245e-17         0.000901
+    33     2.609265e-08     5.290603e-17         0.000939
+    34     2.530975e-08     1.328510e-18         0.000906
 Results of Nonlinear Solver Algorithm
  * Algorithm: Modified Powell
  * Starting Point: [-1.0]
@@ -72,10 +72,10 @@
  * Inf-norm of residuals: 0.000000
  * Convergence: true
  * Message: algorithm estimates that the relative error between x and the solution is at most tol
- * Total time: 0.440380 seconds
+ * Total time: 0.458293 seconds
  * Function Calls: 34
  * Jacobian Calls (df/dx): 2
sol = ϕ((t0, tf), x0, p0_sol.x)                             # get the optimal trajectory
-plot(sol)                                                   # plot
Example block output

Automatic differentiation (wrong way)

Now, we want to provide $J_S$ to the solver, thanks to the $\texttt{ForwardDiff.jl}$ package. This Jacobian is computed with the variational equation, and leads to a false result in our case.

Details. ( To be removed )

Denoting $z = (x,p)$, we have

\[ \varphi(t_0, z_0, t_f) = z_0 + \int_{t_0}^{t_f} \vec H\big(\varphi(t_0, z_0, t)\big) \,\mathrm dt. \]

If we assume that $z_0 \to \varphi(t_0, z_0, t_f)$ is differentiable, we have

\[ \frac{\partial \varphi}{\partial z_0}(t_0, z_0, t_f)\cdot \delta z_0 = \delta z_0 + \int_{t_0}^{t_f} \vec H'\big(\varphi(t_0, z_0, t)\big)\cdot \left( \frac{\partial \varphi}{\partial z_0}(t_0, z_0, t) \cdot \delta z_0 \right) \,\mathrm dt, \]

and so, $z_0 \to \frac{\partial \varphi}{\partial z_0}(t_0, z_0, t_f)\cdot \delta z_0$ is solution of the variational equations

\[ \frac{\partial \delta z}{\partial t}(t) = \vec H'\big(\varphi(t_0, z_0, t_f)\big) \cdot \delta z(t), \qquad \delta z(t_0) = \delta z_0.\]

In the studied optimal control problem, we have

\[ \vec H(x,p) = (\mathrm{sign}(p), -1) \]

and so, we have $\vec H'(z) = 0_2$ almost everywhere, which implies

\[ \frac{\partial \varphi}{\partial z_0}(t_0, z_0, t_f) \cdot \delta z_0 = \mathrm{exp}\big((t_f-t_0) 0_2 \big)\cdot \delta z_0 = \delta z_0.\]

The Jacobian of the shooting function is then given by

\[ S'(p_0) = \pi \left( \frac{\partial \varphi}{\partial p_0}(t_0, x_0, p_0, t_f) \right) = \pi \left( \frac{\partial \varphi}{\partial z_0}(t_0, z_0, t_f) \cdot (0,1) \right) = \pi(0,1) = 0. \]

+plot(sol) # plotExample block output

Automatic differentiation (wrong way)

Now, we want to provide $J_S$ to the solver, thanks to the $\texttt{ForwardDiff.jl}$ package. This Jacobian is computed with the variational equation, and leads to a false result in our case.

Details. ( To be removed )

Denoting $z = (x,p)$, we have

\[ \varphi(t_0, z_0, t_f) = z_0 + \int_{t_0}^{t_f} \vec H\big(\varphi(t_0, z_0, t)\big) \,\mathrm dt. \]

If we assume that $z_0 \to \varphi(t_0, z_0, t_f)$ is differentiable, we have

\[ \frac{\partial \varphi}{\partial z_0}(t_0, z_0, t_f)\cdot \delta z_0 = \delta z_0 + \int_{t_0}^{t_f} \vec H'\big(\varphi(t_0, z_0, t)\big)\cdot \left( \frac{\partial \varphi}{\partial z_0}(t_0, z_0, t) \cdot \delta z_0 \right) \,\mathrm dt, \]

and so, $z_0 \to \frac{\partial \varphi}{\partial z_0}(t_0, z_0, t_f)\cdot \delta z_0$ is solution of the variational equations

\[ \frac{\partial \delta z}{\partial t}(t) = \vec H'\big(\varphi(t_0, z_0, t_f)\big) \cdot \delta z(t), \qquad \delta z(t_0) = \delta z_0.\]

In the studied optimal control problem, we have

\[ \vec H(x,p) = (\mathrm{sign}(p), -1) \]

and so, we have $\vec H'(z) = 0_2$ almost everywhere, which implies

\[ \frac{\partial \varphi}{\partial z_0}(t_0, z_0, t_f) \cdot \delta z_0 = \mathrm{exp}\big((t_f-t_0) 0_2 \big)\cdot \delta z_0 = \delta z_0.\]

The Jacobian of the shooting function is then given by

\[ S'(p_0) = \pi \left( \frac{\partial \varphi}{\partial p_0}(t_0, x_0, p_0, t_f) \right) = \pi \left( \frac{\partial \varphi}{\partial z_0}(t_0, z_0, t_f) \cdot (0,1) \right) = \pi(0,1) = 0. \]

Details. @@ -132,17 +132,17 @@ p0_sol = fsolve(S!, JS!, ξ, show_trace = true) # solve println(p0_sol)
Iter     f(x) inf-norm    Step 2-norm      Step time
 ------   --------------   --------------   --------------
-     1     3.000000e+00     0.000000e+00         0.006296
-     2     5.000000e+00     1.000000e+04         0.059590
-     3     5.000000e+00     3.906250e+03         0.000138
-     4     5.000000e+00     3.906250e+03         0.001273
-     5     5.000000e+00     3.906250e+03         0.000097
-     6     5.000000e+00     5.493164e+02         0.000108
-     7     5.000000e+00     7.724762e+01         0.000110
-     8     9.550781e-01     1.086295e+01         0.001112
-     9     2.633355e-08     2.280435e-01         0.000912
-    10     4.034347e-08     1.733639e-16         0.000973
-    11     7.549755e-09     6.346772e-17         0.051010
+     1     3.000000e+00     0.000000e+00         0.006635
+     2     5.000000e+00     1.000000e+04         0.063100
+     3     5.000000e+00     3.906250e+03         0.000178
+     4     5.000000e+00     3.906250e+03         0.001448
+     5     5.000000e+00     3.906250e+03         0.000107
+     6     5.000000e+00     5.493164e+02         0.000086
+     7     5.000000e+00     7.724762e+01         0.000108
+     8     9.550781e-01     1.086295e+01         0.001187
+     9     2.633355e-08     2.280435e-01         0.000989
+    10     4.034347e-08     1.733639e-16         0.001017
+    11     7.549755e-09     6.346772e-17         0.001007
 Results of Nonlinear Solver Algorithm
  * Algorithm: Modified Powell (User Jac, Expert)
  * Starting Point: [-1.0]
@@ -150,10 +150,10 @@
  * Inf-norm of residuals: 0.000000
  * Convergence: true
  * Message: algorithm estimates that the relative error between x and the solution is at most tol
- * Total time: 0.121820 seconds
+ * Total time: 0.075877 seconds
  * Function Calls: 11
  * Jacobian Calls (df/dx): 2
sol = ϕ((t0, tf), x0, p0_sol.x)                             # get the optimal trajectory
-plt = plot(sol)                                             # plot
Example block output

Automatic differentiation (good way)

The goal is to provide the true Jacobian of $S$ by using the $\texttt{ForwardDiff}$ package, and so we need to indicate to the solver that the dynamic of the system change when $p = 0$.

To understang why we need to give this information to the solver, see the following details.

Details. (To be removed)

The problem is that the Hamiltonian $H$ is not differentiable everywhere due to the maximizing control. This control is bang-bang ($u = 1$ and $u = -1$).

Let now construct the two smooth Hamiltonians associated to these two controls

\[ H^+(x,p) = h(x,p,-1,1) = -x + p \qquad \text{and} \qquad H^-(x,p) = h(x,p,-1,-1) = -x - p. \]

Their associated vector fields are given by

\[ \vec H^+(x,p) = (1,1) \qquad \text{and} \qquad \vec H^-(x,p) = (-1, 1), \]

and their associated flow correspond to

\[ \varphi^+(t_0, z_0, t_f) = z_0 + \left( \begin{array}{c} 1 \\ 1 \end{array} \right) (t_f -t_0) +plt = plot(sol) # plotExample block output

Automatic differentiation (good way)

The goal is to provide the true Jacobian of $S$ by using the $\texttt{ForwardDiff}$ package, and so we need to indicate to the solver that the dynamic of the system change when $p = 0$.

To understang why we need to give this information to the solver, see the following details.

Details. (To be removed)

The problem is that the Hamiltonian $H$ is not differentiable everywhere due to the maximizing control. This control is bang-bang ($u = 1$ and $u = -1$).

Let now construct the two smooth Hamiltonians associated to these two controls

\[ H^+(x,p) = h(x,p,-1,1) = -x + p \qquad \text{and} \qquad H^-(x,p) = h(x,p,-1,-1) = -x - p. \]

Their associated vector fields are given by

\[ \vec H^+(x,p) = (1,1) \qquad \text{and} \qquad \vec H^-(x,p) = (-1, 1), \]

and their associated flow correspond to

\[ \varphi^+(t_0, z_0, t_f) = z_0 + \left( \begin{array}{c} 1 \\ 1 \end{array} \right) (t_f -t_0) \qquad \text{and} \qquad \varphi^-(t_0, z_0, t_f) = z_0 + \left( \begin{array}{c} -1 \\ \phantom{-} 1 \end{array} \right) (t_f -t_0).\]

If we assume that the optimal structure of the problem is negative then positive bangs, then the associated flow is defined by

\[ \varphi(t_0, z_0, t_f) = \varphi^+ \big( t_1(z_0), \varphi^-\big(t_0, z_0, t_1(z_0)\big), t_f \big),\]

with the following condition

\[ \pi_p \big( \varphi^-(t_0, z_0, t_1(z_0)) \big) = 0,\]

where $\pi_p(x,p) = p$ is the classical $p$-space projection. By devlopping this last condition, an explicit form of the function $t_1(\cdot)$ is given by

\[ t_1(x_0, p_0) = t_0 - p_0.\]

Finally, we have

\[ \begin{align*} \frac{\partial \varphi}{\partial z_0} @@ -249,25 +249,25 @@ JS(ξ) : 2.0000000000000018 Iter f(x) inf-norm Step 2-norm Step time ------ -------------- -------------- -------------- - 1 3.000000e+00 0.000000e+00 1.318328 - 2 6.994405e-15 2.250000e+00 0.007805 - 3 4.551914e-15 1.262177e-29 0.000191 + 1 3.000000e+00 0.000000e+00 1.394785 + 2 6.994405e-15 2.250000e+00 0.008386 + 3 4.551914e-15 1.262177e-29 0.000203 4 3.441691e-15 1.774937e-30 0.000142 - 5 2.775558e-15 1.972152e-31 0.000141 - 6 3.441691e-15 1.972152e-31 0.000153 - 7 2.775558e-15 1.972152e-31 0.000141 - 8 2.997602e-15 1.774937e-30 0.000294 - 9 1.443290e-15 1.972152e-31 0.000141 - 10 2.997602e-15 7.888609e-31 0.000128 - 11 2.997602e-15 1.972152e-31 0.000135 - 12 2.997602e-15 1.972152e-31 0.000276 - 13 2.997602e-15 1.972152e-31 0.000143 - 14 1.443290e-15 1.972152e-31 0.000128 - 15 2.343750e-02 1.373291e-04 0.000135 - 16 2.997602e-15 1.373291e-04 0.000129 - 17 2.997602e-15 1.972152e-31 0.000143 - 18 1.443290e-15 1.972152e-31 0.000128 - 19 1.464844e-03 5.364418e-07 0.000137 + 5 2.775558e-15 1.972152e-31 0.000148 + 6 3.441691e-15 1.972152e-31 0.000133 + 7 2.775558e-15 1.972152e-31 0.000140 + 8 2.997602e-15 1.774937e-30 0.000300 + 9 1.443290e-15 1.972152e-31 0.000145 + 10 2.997602e-15 7.888609e-31 0.000129 + 11 2.997602e-15 1.972152e-31 0.000139 + 12 2.997602e-15 1.972152e-31 0.000270 + 13 2.997602e-15 1.972152e-31 0.000146 + 14 1.443290e-15 1.972152e-31 0.000129 + 15 2.343750e-02 1.373291e-04 0.000140 + 16 2.997602e-15 1.373291e-04 0.000128 + 17 2.997602e-15 1.972152e-31 0.000141 + 18 1.443290e-15 1.972152e-31 0.000125 + 19 1.464844e-03 5.364418e-07 0.000139 Results of Nonlinear Solver Algorithm * Algorithm: Modified Powell (User Jac, Expert) * Starting Point: [-1.0] @@ -275,7 +275,7 @@ * Inf-norm of residuals: 0.000000 * Convergence: true * Message: iteration is not making good progress, measured by improvement from last 10 iterations - * Total time: 1.328833 seconds + * Total time: 1.405883 seconds * Function Calls: 19 * Jacobian Calls (df/dx): 3

# get optimal trajectory
 sol_ = φ((t0, tf), x0, p0_sol.x[1], saveat=range(t0, tf, 500))
@@ -292,4 +292,4 @@
 plt_u = plot(t, u, label = "u")
 
 plt_xp = plot(plt_x, plt_p, layout=(1, 2))
-plot(plt_xp, plt_u, layout = (2, 1))
Example block output
+plot(plt_xp, plt_u, layout = (2, 1))Example block output
diff --git a/dev/2D-preconditioner-246e2978.svg b/dev/2D-preconditioner-6038766b.svg similarity index 64% rename from dev/2D-preconditioner-246e2978.svg rename to dev/2D-preconditioner-6038766b.svg index 38d4e4a..54b394d 100644 --- a/dev/2D-preconditioner-246e2978.svg +++ b/dev/2D-preconditioner-6038766b.svg @@ -1,2107 +1,2107 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/2D-preconditioner-2ee2be80.svg b/dev/2D-preconditioner-a1c3d910.svg similarity index 86% rename from dev/2D-preconditioner-2ee2be80.svg rename to dev/2D-preconditioner-a1c3d910.svg index 808b9f6..57cc261 100644 --- a/dev/2D-preconditioner-2ee2be80.svg +++ b/dev/2D-preconditioner-a1c3d910.svg @@ -1,194 +1,194 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/2D-preconditioner-e66ead1f.svg b/dev/2D-preconditioner-d0a533dd.svg similarity index 85% rename from dev/2D-preconditioner-e66ead1f.svg rename to dev/2D-preconditioner-d0a533dd.svg index 01ba0df..a757940 100644 --- a/dev/2D-preconditioner-e66ead1f.svg +++ b/dev/2D-preconditioner-d0a533dd.svg @@ -1,190 +1,190 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/2D-preconditioner-b8e7ebd5.svg b/dev/2D-preconditioner-d454d92a.svg similarity index 81% rename from dev/2D-preconditioner-b8e7ebd5.svg rename to dev/2D-preconditioner-d454d92a.svg index df99e8e..a508536 100644 --- a/dev/2D-preconditioner-b8e7ebd5.svg +++ b/dev/2D-preconditioner-d454d92a.svg @@ -1,146 +1,146 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/2D-preconditioner.html b/dev/2D-preconditioner.html index f0b9d97..eb908db 100644 --- a/dev/2D-preconditioner.html +++ b/dev/2D-preconditioner.html @@ -93,23 +93,23 @@ p0_sol = fsolve(S₂!, JS₂!, ξ, show_trace = true) # solve println(p0_sol) # print solution
Iter     f(x) inf-norm    Step 2-norm      Step time
 ------   --------------   --------------   --------------
-     1     3.845299e+00     0.000000e+00         0.003356
-     2     5.000000e+00     1.559496e+00         4.607092
-     3     5.000000e+00     4.983094e-01         0.000139
-     4     5.000000e+00     4.983094e-01         0.000371
-     5     5.000000e+00     4.983094e-01         0.000088
-     6     2.825589e+00     9.417471e-02         0.000252
-     7     3.275582e+00     5.569993e-02         0.000262
-     8     1.835694e+00     1.605481e-02         0.000189
-     9     5.000000e+00     4.108358e-02         0.000080
-    10     8.767661e-01     2.198078e-02         0.000172
-    11     1.052824e+00     2.476848e-03         0.000242
-    12     2.252319e-01     7.373618e-04         0.000174
-    13     6.733089e-02     6.111140e-05         0.000268
-    14     4.061534e-03     3.236777e-06         0.000173
-    15     6.974839e-05     1.047584e-08         0.000240
-    16     7.333033e-08     3.198326e-12         0.000172
-    17     1.322831e-12     3.527840e-18         0.000244
+     1     3.845299e+00     0.000000e+00         0.003597
+     2     5.000000e+00     1.559496e+00         5.080205
+     3     5.000000e+00     4.983094e-01         0.000082
+     4     5.000000e+00     4.983094e-01         0.000382
+     5     5.000000e+00     4.983094e-01         0.000028
+     6     2.825589e+00     9.417471e-02         0.000239
+     7     3.275582e+00     5.569993e-02         0.000178
+     8     1.835694e+00     1.605481e-02         0.000176
+     9     5.000000e+00     4.108358e-02         0.000023
+    10     8.767661e-01     2.198078e-02         0.000157
+    11     1.052824e+00     2.476848e-03         0.000165
+    12     2.252319e-01     7.373618e-04         0.000157
+    13     6.733089e-02     6.111140e-05         0.000164
+    14     4.061534e-03     3.236777e-06         0.000158
+    15     6.974839e-05     1.047584e-08         0.000169
+    16     7.333033e-08     3.198326e-12         0.000157
+    17     1.322831e-12     3.527840e-18         0.000167
 Results of Nonlinear Solver Algorithm
  * Algorithm: Modified Powell (User Jac, Expert)
  * Starting Point: [-0.5]
@@ -117,7 +117,7 @@
  * Inf-norm of residuals: 0.000000
  * Convergence: true
  * Message: algorithm estimates that the relative error between x and the solution is at most tol
- * Total time: 4.613529 seconds
+ * Total time: 5.086217 seconds
  * Function Calls: 17
  * Jacobian Calls (df/dx): 2
sol = ϕ((t0, tf), x0, [η(p0_sol.x[1]), p0_sol.x[1]],
     saveat=range(t0, tf, 500))                                  # get optimal trajectory
@@ -137,7 +137,7 @@
 plt_u  = plot(t, u,  label = "u" )
 
 plt_xp = plot(plt_x⁰, plt_p⁰, plt_x, plt_p, layout=(2, 2))
-plot(plt_xp, plt_u, layout = grid(2,1, heights = [2/3, 1/3]), size=(700, 600))
Example block output

Construction of the geometric preconditioner

The goal is now to use the geometric preconditioning method proposed in [mettre article]. For this purpose, the first step is to create some points in the boundary of the accessible augmented set, and to fit an ellipse on these points.

The second step is to create the linear diffeomorphism $\phi \colon \mathbb R^2 \to \mathbb R^2, (x^0,x_1) \to Ax + B$ that transforms the fitted ellipse into the unit circle, and that satisfy the condition

\[ \frac{\partial \phi}{\partial x^0} = k e_1, \]

with $k>0$ and $e_1 = (1,0)$. In this context, we denote

\[ A = \left( \begin{array}{cc} k & A_{x^0} \\ 0 & A_{x_1} \end{array} \right) +plot(plt_xp, plt_u, layout = grid(2,1, heights = [2/3, 1/3]), size=(700, 600))Example block output

Construction of the geometric preconditioner

The goal is now to use the geometric preconditioning method proposed in [mettre article]. For this purpose, the first step is to create some points in the boundary of the accessible augmented set, and to fit an ellipse on these points.

The second step is to create the linear diffeomorphism $\phi \colon \mathbb R^2 \to \mathbb R^2, (x^0,x_1) \to Ax + B$ that transforms the fitted ellipse into the unit circle, and that satisfy the condition

\[ \frac{\partial \phi}{\partial x^0} = k e_1, \]

with $k>0$ and $e_1 = (1,0)$. In this context, we denote

\[ A = \left( \begin{array}{cc} k & A_{x^0} \\ 0 & A_{x_1} \end{array} \right) \qquad \text{and} \qquad B = \left( \begin{array}{c} B_{x^0} \\ B_{x_1} \end{array} \right).\]

This diffeomorphism is given from the semi-axis $a,b >0$, the angle $\theta \in [0, \frac{\pi}{2}[$ between the semi-axis $b$ and the $x$-axis, and the center $c \in \mathbb R^2$ by

\[ \phi(x) = r(-\beta_0) s(a^{-1}, b^{-1}) r(\theta) (x - c),\]

where $r$ and $s$ correspond respectively to the rotation and the scale matrix, defined by

\[ r(\theta) = \left( \begin{array}{cc} \phantom - \cos(\theta) & \sin(\theta) \\ -\sin(\theta) & \cos(\theta) \end{array} \right) \qquad \text{and} \qquad @@ -202,7 +202,7 @@ plot!(yₑ[2,:], yₑ[1,:], label = "") plot!( xlabel = "y", ylabel = "y⁰") -plot(plt_x, plt_y, layout = (1,2), size=(800, 400))Example block output

The general shooting function $T \colon \mathbb R^2 \to \mathbb R$ in the new coordinate is defined by

\[ T(q) = A_{x_1} \varphi \big(t_0, x_0, p_0(A^\top q), t_f) + B_{x_1} - y_T,\]

where the function $p_0 \colon \mathbb R^2 \to \mathbb R^2$ correspond to the mapping between the final and the initial augmented costate, and is given by

\[ p_0(p^0, p_1) = (p^0, p_1 + p^0 t_f),\]

and $y_T = A_{x_1} x_T + B_x$ to the target in the new system of coordinates. By using the definition of $S$ and $y_T$, we obtain

\[\begin{align*} +plot(plt_x, plt_y, layout = (1,2), size=(800, 400))Example block output

The general shooting function $T \colon \mathbb R^2 \to \mathbb R$ in the new coordinate is defined by

\[ T(q) = A_{x_1} \varphi \big(t_0, x_0, p_0(A^\top q), t_f) + B_{x_1} - y_T,\]

where the function $p_0 \colon \mathbb R^2 \to \mathbb R^2$ correspond to the mapping between the final and the initial augmented costate, and is given by

\[ p_0(p^0, p_1) = (p^0, p_1 + p^0 t_f),\]

and $y_T = A_{x_1} x_T + B_x$ to the target in the new system of coordinates. By using the definition of $S$ and $y_T$, we obtain

\[\begin{align*} T(q) &= A_{x_1} \varphi \big(t_0, x_0, p_0(A^\top q), t_f) + B_{x_1} - (A_{x_1} x_T + B_{x_1}) \\ &= A_{x_1} (S \circ p_0)(A^\top q). \end{align*}\]

which highlight that the proposed geometric preconditioning method is a left and right side preconditioner of the shooting function. Finally, we define the two shooting functions $T_1$ and $T_2$ by using the two method of normalization used before for the function $S$

\[ T_1(p) = T(-1, p) \qquad \text{and} \qquad T_2(p) = T \big(\eta(p), p \big).\]

p₀(p) = [p[1], p[2] + p[1]*tf]                                  # function p₀(⋅)
@@ -248,24 +248,24 @@
 q_sol = fsolve(T₂!, JT₂!, ξ, show_trace = true)                 # solve
 println(q_sol)                                                  # print solution
Iter     f(x) inf-norm    Step 2-norm      Step time
 ------   --------------   --------------   --------------
-     1     3.223430e-01     0.000000e+00         0.000415
-     2     7.034101e-02     1.406250e-01         4.899762
-     3     1.135104e-02     1.095650e-02         0.000259
-     4     1.033639e-04     4.056828e-04         0.000186
-     5     2.155036e-08     3.426091e-08         0.000179
-     6     8.398089e-16     1.489880e-15         0.000168
-     7     3.119290e-16     2.262576e-30         0.000175
-     8     2.399454e-17     1.659615e-31         0.000167
-     9     2.399454e-17     1.152511e-33         0.000176
-    10     3.119290e-16     1.947743e-31         0.000165
-    11     2.399454e-17     2.028246e-31         0.000458
-    12     2.399454e-17     2.097361e-32         0.000177
-    13     2.399454e-17     2.593149e-33         0.000173
-    14     2.399454e-17     6.482873e-34         0.000165
-    15     2.399454e-17     1.620718e-34         0.000170
-    16     2.399454e-17     4.051796e-35         0.000185
-    17     2.399454e-17     1.012949e-35         0.000174
-    18     2.399454e-17     2.532372e-36         0.000164
+     1     3.223430e-01     0.000000e+00         0.000436
+     2     7.034101e-02     1.406250e-01         5.067886
+     3     1.135104e-02     1.095650e-02         0.000633
+     4     1.033639e-04     4.056828e-04         0.000271
+     5     2.155036e-08     3.426091e-08         0.000222
+     6     8.398089e-16     1.489880e-15         0.000206
+     7     3.119290e-16     2.262576e-30         0.000178
+     8     2.399454e-17     1.659615e-31         0.000160
+     9     2.399454e-17     1.152511e-33         0.000167
+    10     3.119290e-16     1.947743e-31         0.000160
+    11     2.399454e-17     2.028246e-31         0.000449
+    12     2.399454e-17     2.097361e-32         0.000166
+    13     2.399454e-17     2.593149e-33         0.000167
+    14     2.399454e-17     6.482873e-34         0.000156
+    15     2.399454e-17     1.620718e-34         0.000162
+    16     2.399454e-17     4.051796e-35         0.000156
+    17     2.399454e-17     1.012949e-35         0.000164
+    18     2.399454e-17     2.532372e-36         0.000157
 Results of Nonlinear Solver Algorithm
  * Algorithm: Modified Powell (User Jac, Expert)
  * Starting Point: [0.5]
@@ -273,7 +273,7 @@
  * Inf-norm of residuals: 0.000000
  * Convergence: true
  * Message: iteration is not making good progress, measured by improvement from last 10 iterations
- * Total time: 4.903338 seconds
+ * Total time: 5.071917 seconds
  * Function Calls: 18
  * Jacobian Calls (df/dx): 2
p0_sol = p₀(transpose(A)*[η(q_sol.x[1]), q_sol.x[1]])           # get the optimal initial costate in old coordinates
 sol = ϕ((t0, tf), x0, p0_sol, saveat=range(t0, tf, 500))        # get optimal trajectory
@@ -293,7 +293,7 @@
 plt_u  = plot(t, u,  label = "u" )
 
 plt_xp = plot(plt_x⁰, plt_p⁰, plt_x, plt_p, layout=(2, 2))
-plot(plt_xp, plt_u, layout = grid(2,1, heights = [2/3, 1/3]), size=(700, 600))
Example block output

Comparison

It is shown in [mettre article] that if the boundary of augmented accessible set is the fitted ellipse then the shooting function $T_2$ is the identity function. Due to the error of the approximation, the function $T_2$ is not the identity, but we hope that this function is close to this ideal function, and so that the convergence of $T_2$ is faster than the one of $S_2$.

The following code compares the convergence of these two shooting function. We also study for these functions the domain of initial guesses for the solver that hit the bounds $[-1, 1]$ during the solving process. 

## initial guesses
+plot(plt_xp, plt_u, layout = grid(2,1, heights = [2/3, 1/3]), size=(700, 600))
Example block output

Comparison

It is shown in [mettre article] that if the boundary of augmented accessible set is the fitted ellipse then the shooting function $T_2$ is the identity function. Due to the error of the approximation, the function $T_2$ is not the identity, but we hope that this function is close to this ideal function, and so that the convergence of $T_2$ is faster than the one of $S_2$.

The following code compares the convergence of these two shooting function. We also study for these functions the domain of initial guesses for the solver that hit the bounds $[-1, 1]$ during the solving process. 

## initial guesses
 N = 1000; ξ = range(-1.,1.,N)
 
 # initialization of the matrix
@@ -368,4 +368,4 @@
 scatter!(1,1, markercolor = :red, markerstrokecolor = :red, marker = 2, label = "not converged")
 
 plt_2 = plot(plt_21, plt_22, layout = (1,2))
-plot(plt_1, plt_2, layout = grid(2,1, heights = [0.5, 0.5]), size=(700, 600))
Example block output +plot(plt_1, plt_2, layout = grid(2,1, heights = [0.5, 0.5]), size=(700, 600))Example block output diff --git a/dev/index.html b/dev/index.html index 1e1dc6f..3a4805d 100644 --- a/dev/index.html +++ b/dev/index.html @@ -4,7 +4,7 @@ \text{s.c.}~\dot x(t) = u(t), & t\in [t_0, t_f]~\mathrm{a.e.}, \\[0.5em] \phantom{\mathrm{s.c.}~} u(t) \in [-1,1], & t\in [t_0, t_f], \\[0.5em] \phantom{\mathrm{s.c.}~} x(t_0) = x_0, \quad x(t_f) = x_f, - \end{array} \right.\]

with $x_0$, $t_0$, $x_f$ and $t_f$ fixed. This problem is simple, and can be analytically solve without the use of numerical method. However, the goal is to solve this problem by indirect shooting.

Dependencies

All the numerical simulations to generate this documentation from MRI.jl are performed with the following packages.

using Pkg
+    \end{array} \right.\]
with $x_0$, $t_0$, $x_f$ and $t_f$ fixed. This problem is simple, and can be analytically solve without the use of numerical method. However, the goal is to solve this problem by indirect shooting.  
Dependencies
All the numerical simulations to generate this documentation are performed with  the following packages.
using Pkg
 Pkg.status()
Status `~/work/preconditioning/preconditioning/docs/Project.toml`
   [0c46a032] DifferentialEquations v7.13.0
   [e30172f5] Documenter v1.5.0
@@ -12,4 +12,4 @@
   [4fb4c77a] GeometricPreconditioner v0.1.0 `~/work/preconditioning/preconditioning`
   [4854310b] MINPACK v1.3.0
   [5f98b655] OptimalControl v0.11.0
-  [91a5bcdd] Plots v1.40.5

+  [91a5bcdd] Plots v1.40.5
diff --git a/dev/search_index.js b/dev/search_index.js index 79e7029..219a5e8 100644 --- a/dev/search_index.js +++ b/dev/search_index.js @@ -1,3 +1,3 @@ var documenterSearchIndex = {"docs": -[{"location":"2D-preconditioner.html#Augmented-formulation","page":"Geometric preconditioner","title":"Augmented formulation","text":"","category":"section"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"We consider the augmented formulation of the same problem as before","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":" left beginarrayll\n displaystyle min_xu x^0(t_f) 1em\n textscdot x^0(t) = x_1(t) tin t_0 t_fmathrmae 05em\n phantommathrmsc dot x_1(t) = u(t) tin t_0 t_fmathrmae 05em\n phantommathrmsc u(t) in -11 tin t_0 t_f 05em\n phantommathrmsc x(t_0) = x_0 quad x_1(t_f) = x_f\n endarray right","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"with fixed x_0, t_0, x_f and t_f, and where x = (x^0 x_1). The augmented Hamiltonian of this problem is given by ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":" H(xp) = p^0 x_1 + lvert p_1 rvert","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"where p = (p^0 p_1) is the augmented costate, composed by the costate p_1 associated to the state x_1 and multiplier p^0 associated to the cost x^0 which is no longer fixed to -1. The flow varphi associated to the true Hamiltonian H is computed thanks to the textttflow and the textttHamiltonian functions.","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"The general shooting function S colon mathbb R^2 to mathbb R associated to this Hamiltonian is defined by ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":" S(p_0) = pi_x big( varphi(t_0 x_0 p_0 t_f) big)","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"where pi_x(x^0 x_1 p^0 x_1) = x_1 is still the state projection. ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"We are interested in two normalization of this shooting function, thanks to the homogeneity of the BC-extremals on the augmented costate, which are ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":" S_1(p_0) = S(-1 p_0) qquad textand qquad S_2(p_0) = S big( eta(p) p_0 big)","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"where the function eta colon -1 1 to mathbb R is defined by ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":" eta(p) = -sqrt 1 - p^2","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"The following code creates and plots these shooting functions. ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"using OptimalControl\nusing Plots\nusing ForwardDiff\nusing DifferentialEquations\nusing MINPACK\nusing Statistics\n\nglobal α = 1\n\nfunction condition(z, t, integrator) # event when condition(z,t,integrator) == 0\n x⁰,x,p⁰,p = z\n return p\nend\n\nfunction affect!(integrator) # action when condition == 0 \n global α = -α\n nothing\nend\n\ncb = ContinuousCallback(condition, affect!) # callback \nH(x, p) = p[1] * x[2] + α * p[2]\nϕ_ = Flow(Hamiltonian(H), callback = cb) # flow with maximizing control \n\nfunction ϕ(t0, x0, p0, tf; kwargs...) \n if p0[2] == 0\n global α = -sign(p0[1])\n else\n global α = sign(p0[2])\n end\n return ϕ_(t0, x0, p0, tf; kwargs...)\nend\n\nfunction ϕ((t0, tf), x0, p0; kwargs...) # flow for plot\n if p0[2] == 0\n global α = -sign(p0[1])\n else\n global α = sign(p0[2])\n end\n return ϕ_((t0, tf), x0, p0; kwargs...)\nend\n\nt0 = 0 # initial time\nx0 = [0,0] # initial augmented state\ntf = 5 # final time\nxT = 0 # final state\n\nπ((x,p)) = x[2] # projection on state space\nη(p0) = -sqrt.(1 - p0.^2) # function η(⋅)\n\nS(p0) = π( ϕ(t0, x0, p0, tf) ) - xT # general hooting function\n\nS₁(p0) = S([-1, p0]) # normalization 1\nS₂(p0) = abs(p0) < 1 ? S([η(p0), p0]) : sign(p0)*tf - xT # normalization 2\n\n# Plot\nplt_S1 = plot(range(-7, 2, 500), S₁ , color = 2, label = \"\")\nplot!(xlabel = \"p₀\", ylabel = \"S₁\", xlim = [-7,2])\n\nplt_S2 = plot(range(-1, 1, 500), S₂, color = 3)\nplot!(xlabel = \"p₀\", ylabel = \"S₂\", xlim = [-1,1], legend=false)\n\nS_(p⁰, p) = S([p⁰, p])\nplt_S = surface(range(0, -1.5, 100), range(-8, 2, 100), S_, camera = (30,30))\nsurface!(xlabel = \"p⁰\", ylabel = \"p₀\", zlabel = \"φₓ\", xflip = true)\nplot3d!(-1*ones(100), range(-8, 2, 100), S₁.(range(-8, 2, 100)), label = \"S₁\")\nplot3d!(η.(range(-1, 1, 100)), range(-1, 1, 100), S₂.(range(-1, 1, 100)), label = \"S₂\")\n\nplt_S12 = plot(plt_S1, plt_S2, layout = (1,2))\nplt_total = plot(plt_S, plt_S12, layout = grid(2,1, heights = [2/3, 1/3]), size=(700, 600))\n\n@gif for i ∈ [range(30, 90, 50); 90*ones(25); range(90, 30, 50); 30*ones(25)]\n plot!(plt_total[1], camera=(i,i), \n zticks = i==90 ? false : true,\n zlabel = i==90 ? \"\" : \"S\" )\nend","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"As done before, we use the solver texttthybrd1 from the textttMINPACKjl package to find a zero of S_2.","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"global iterate_S2 = Vector{Float64}() # global vector to store iterates of the solver\n\nfunction S₂!(s₂, ξ) # intermediate function\n push!(iterate_S2, ξ[1]) # stock the iterate\n return (s₂[:] .= S₂(ξ[1]); nothing) \nend\n\nJS₂(ξ) = ForwardDiff.jacobian(p0 -> [S₂(p0[1])], ξ) # compute jacobian by forward differentiation\nJS₂!(js₂, ξ) = (js₂[:] .= JS₂(ξ); nothing) # intermediate function\n\nξ = [-0.5] # initial guess\np0_sol = fsolve(S₂!, JS₂!, ξ, show_trace = true) # solve\nprintln(p0_sol) # print solution","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"sol = ϕ((t0, tf), x0, [η(p0_sol.x[1]), p0_sol.x[1]], \n saveat=range(t0, tf, 500)) # get optimal trajectory\n\n# plot\nt = sol.t\nx⁰ = [sol.u[i][1] for i in 1:length(sol.u)]\nx = [sol.u[i][2] for i in 1:length(sol.u)]\np⁰ = [sol.u[i][3] for i in 1:length(sol.u)]\np = [sol.u[i][4] for i in 1:length(sol.u)]\nu = sign.(p)\n\nplt_x⁰ = plot(t, x⁰, label = \"x⁰\")\nplt_x = plot(t, x , label = \"x\" )\nplt_p⁰ = plot(t, p⁰, label = \"p⁰\", ylim=[-1,1])\nplt_p = plot(t, p , label = \"p\" )\nplt_u = plot(t, u, label = \"u\" )\n\nplt_xp = plot(plt_x⁰, plt_p⁰, plt_x, plt_p, layout=(2, 2))\nplot(plt_xp, plt_u, layout = grid(2,1, heights = [2/3, 1/3]), size=(700, 600))","category":"page"},{"location":"2D-preconditioner.html#Construction-of-the-geometric-preconditioner","page":"Geometric preconditioner","title":"Construction of the geometric preconditioner","text":"","category":"section"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"The goal is now to use the geometric preconditioning method proposed in [mettre article]. For this purpose, the first step is to create some points in the boundary of the accessible augmented set, and to fit an ellipse on these points. ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"The second step is to create the linear diffeomorphism phi colon mathbb R^2 to mathbb R^2 (x^0x_1) to Ax + B that transforms the fitted ellipse into the unit circle, and that satisfy the condition ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":" fracpartial phipartial x^0 = k e_1 ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"with k0 and e_1 = (10). In this context, we denote ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":" A = left( beginarraycc k A_x^0 0 A_x_1 endarray right)\n qquad textand qquad \n B = left( beginarrayc B_x^0 B_x_1 endarray right)","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"This diffeomorphism is given from the semi-axis ab 0, the angle theta in 0 fracpi2 between the semi-axis b and the x-axis, and the center c in mathbb R^2 by ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":" phi(x) = r(-beta_0) s(a^-1 b^-1) r(theta) (x - c)","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"where r and s correspond respectively to the rotation and the scale matrix, defined by ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":" r(theta) = left( beginarraycc phantom - cos(theta) sin(theta) -sin(theta) cos(theta) endarray right)\n qquad textand qquad\n s(ab) = left( beginarraycc a 0 0 b endarray right)","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"and where beta_0 = arctan left(fraca sin(theta)b cos(theta) right). ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"
\n
\n \n fit_ellipse — Function\n
\n
function fit_ellipse(x, y)\n    M = hcat(x.^2, x.*y, y.^2, x, y)                            # quadratic form of ellipse\n    p = M\\ones(length(x))                                       # fit parameters for the ellipse\n    A, B, C, D, E = p\n    F = -1.0\n    # calculate the parameters from quadratic\n    Δ = B^2 - 4*A*C\n    Λ = (A-C)^2 + B^2\n    b, a = [-sqrt(clamp( 2*(A*E^2 + C*D^2 - B*D*E + Δ*F)*\n            ( (A+C) + op(sqrt(Λ)) ), 0, Inf)) / Δ   for op in (+, -)]\n    θ = atan(-B, C-A)/2\n    c = [(2*C*D - B*E)/Δ, (2*A*E - B*D)/Δ]\n    return a, b, -θ+Base.π/2, c\nend
\n
\n
","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"function fit_ellipse(x, y) # hide\n M = hcat(x.^2, x.*y, y.^2, x, y) # hide\n p = M\\ones(length(x)) # hide\n A, B, C, D, E = p # hide\n F = -1.0 # hide\n Δ = B^2 - 4*A*C # hide\n Λ = (A-C)^2 + B^2 # hide\n b, a = [-sqrt(clamp( 2*(A*E^2 + C*D^2 - B*D*E + Δ*F) * ( (A+C) + op(sqrt(Λ)) ), 0, Inf)) / Δ for op in (+, -)] # hide\n θ = atan(-B, C-A)/2 # hide\n c = [(2*C*D - B*E)/Δ, (2*A*E - B*D)/Δ] # hide\n return a, b, -θ+Base.π/2, c # hide\nend # hide\nnothing # hide","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"n = 15 # number of points for fit : 2n\nn_ = 100 # number of points for plot: 2n_\np0 = [[[-1, i] for i ∈ range(-tf, 0, n)]; # initial costate for fit\n [[1, i] for i ∈ range(tf, 0, n)]] \np0_ = [[[-1, i] for i ∈ range(-tf, 0, n_)]; # initial costate for plot\n [[1, i] for i ∈ range(tf, 0, n_)]]\n\nx = zeros(2, 2*n); p = zeros(2,2*n) # init final state and costate\nx_ = zeros(2, 2*n_); p_ = zeros(2, 2*n_)\nfor i = 1:length(p0)\n x[:,i], p[:,i] = ϕ(t0, x0, p0[i], tf) # compute flow for fit\nend\nfor i = 1:length(p0_)\n x_[:,i], p_[:,i] = ϕ(t0, x0, p0_[i], tf) # compute flow for plot\nend\n\na, b, θ, c = fit_ellipse(x[1,:], x[2,:]) # fit ellipse\nr(β) = [[cos(β), sin(β)] [-sin(β), cos(β)]] # 2x2 rotation matrix\ns(a,b) = [[a,0] [0,b]] # 2x2 scale matrix\nβ = range(-Base.π, Base.π; length = 100) # angle for plot ellipse\nxₑ = r(-θ)*s(a,b)*\n transpose(reduce(hcat,[sin.(β), cos.(β)])).+c # points of the ellipse\n\n# construction of the linear diffeomorphism φ(x) = Ax + B\nd = (a*sin(θ))/(b*cos(θ)); β₀ = atan(d) # intermediate values\nA = r(-β₀)*s(1/a,1/b)*r(θ); B = -A*c # calculate A and B\nφ(x) = A*x .+ B # function φ\n\ny = φ(x); y_ = φ(x_); yₑ = φ(xₑ) # compute φ on x, x_ and xₑ\n\n# plot\nplt_x = plot(x_[2,:], x_[1,:], label = \"\")\nscatter!(x[2,:], x[1,:], label=\"Observations\", legend = :topleft)\nplot!(xₑ[2,:], xₑ[1,:], label = \"Fitted ellipse\")\nplot!(xlim = [-15,15], ylim = [-15,15], xlabel = \"x\", ylabel = \"x⁰\")\n\nplt_y = plot(y_[2,:], y_[1,:], label = \"\")\nscatter!(y[2,:], y[1,:], label=\"\")\nplot!(yₑ[2,:], yₑ[1,:], label = \"\")\nplot!( xlabel = \"y\", ylabel = \"y⁰\")\n\nplot(plt_x, plt_y, layout = (1,2), size=(800, 400))","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"The general shooting function T colon mathbb R^2 to mathbb R in the new coordinate is defined by ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":" T(q) = A_x_1 varphi big(t_0 x_0 p_0(A^top q) t_f) + B_x_1 - y_T","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"where the function p_0 colon mathbb R^2 to mathbb R^2 correspond to the mapping between the final and the initial augmented costate, and is given by","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":" p_0(p^0 p_1) = (p^0 p_1 + p^0 t_f)","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"and y_T = A_x_1 x_T + B_x to the target in the new system of coordinates. By using the definition of S and y_T, we obtain ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"beginalign*\n T(q) = A_x_1 varphi big(t_0 x_0 p_0(A^top q) t_f) + B_x_1 - (A_x_1 x_T + B_x_1) \n = A_x_1 (S circ p_0)(A^top q)\nendalign*","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"which highlight that the proposed geometric preconditioning method is a left and right side preconditioner of the shooting function. Finally, we define the two shooting functions T_1 and T_2 by using the two method of normalization used before for the function S ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":" T_1(p) = T(-1 p) qquad textand qquad T_2(p) = T big(eta(p) p big)","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"p₀(p) = [p[1], p[2] + p[1]*tf] # function p₀(⋅)\nAₓ = A[2,2]; Bₓ = B[2] # Aₓ and Bₓ\nyT = Aₓ*xT + Bₓ # target in the new coordinate \nT(q) = Aₓ*(S ∘ p₀)(transpose(A)*q) # state flow in the new coordinates\nT₁(q) = T([-1, q]) # normalization 1\nT₂(q) = abs(q) < 1 ? T([η(q), q]) : sign(q)*(Aₓ*tf + Bₓ) # normalization 2\n\n# Plot\nplt_S1 = plot(range(-3, 3, 500), T₁ , color = 2)\nplot!([-3,3], [yT,yT], color = :black)\nplot!(xlabel = \"q₀\", ylabel = \"T₁\", xlim = [-3,3], legend=false)\n\nplt_S2 = plot(range(-1, 1, 500), T₂, color = 3)\nplot!([-1,1], [yT,yT], color = :black)\nplot!(xlabel = \"q₀\", ylabel = \"T₂\", xlim = [-1,1], legend=false)\n\nT_(q⁰, q) = T([q⁰, q])\nplt_S = surface(range(0, -1.5, 100), range(-3, 3, 100), T_, camera = (30,30))\nsurface!(xlabel = \"q⁰\", ylabel = \"q₀\", zlabel = \"T\", xflip = true)\nplot3d!(-1*ones(100), range(-3, 3, 100), T₁.(range(-3, 3, 100)), label = \"T₁\")\nplot3d!(η.(range(-1, 1, 100)), range(-1, 1, 100), T₂.(range(-1, 1, 100)), label = \"T₂\")\n\nplt_S12 = plot(plt_S1, plt_S2, layout = (1,2))\nplt_total = plot(plt_S, plt_S12, layout = grid(2,1, heights = [2/3, 1/3]), size=(700, 600))\n\n@gif for i ∈ [range(30, 90, 50); 90*ones(25); range(90, 30, 50); 30*ones(25)]\n plot!(plt_total[1], camera=(i,i), \n zticks = i==90 ? false : true,\n zlabel = i==90 ? \"\" : \"T\" )\nend","category":"page"},{"location":"2D-preconditioner.html#Preconditioned-shoot","page":"Geometric preconditioner","title":"Preconditioned shoot","text":"","category":"section"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"global iterate_T2 = Vector{Float64}() # global vector to store iterates of the solver\n\nfunction T₂!(t₂, ξ) # intermediate function\n push!(iterate_T2, ξ[1]) \n return (t₂[:] .= T₂(ξ[1]); nothing) \nend\n\nJT₂(ξ) = ForwardDiff.jacobian(q0 -> [T₂(q0[1])], ξ) # compute jacobian by forward differentiation\nJT₂!(jt₂, ξ) = (jt₂[:] .= JT₂(ξ); nothing) # intermediate function\n\nξ = [0.5] # initial guess\nq_sol = fsolve(T₂!, JT₂!, ξ, show_trace = true) # solve\nprintln(q_sol) # print solution","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"p0_sol = p₀(transpose(A)*[η(q_sol.x[1]), q_sol.x[1]]) # get the optimal initial costate in old coordinates\nsol = ϕ((t0, tf), x0, p0_sol, saveat=range(t0, tf, 500)) # get optimal trajectory\n\n# plot\nt = sol.t\nx⁰ = [sol.u[i][1] for i in 1:length(sol.u)]\nx = [sol.u[i][2] for i in 1:length(sol.u)]\np⁰ = [sol.u[i][3] for i in 1:length(sol.u)]\np = [sol.u[i][4] for i in 1:length(sol.u)]\nu = sign.(p)\n\nplt_x⁰ = plot(t, x⁰, label = \"x⁰\")\nplt_x = plot(t, x , label = \"x\" )\nplt_p⁰ = plot(t, p⁰, label = \"p⁰\", ylim=[-1,1])\nplt_p = plot(t, p , label = \"p\" )\nplt_u = plot(t, u, label = \"u\" )\n\nplt_xp = plot(plt_x⁰, plt_p⁰, plt_x, plt_p, layout=(2, 2))\nplot(plt_xp, plt_u, layout = grid(2,1, heights = [2/3, 1/3]), size=(700, 600))","category":"page"},{"location":"2D-preconditioner.html#Comparison","page":"Geometric preconditioner","title":"Comparison","text":"","category":"section"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"It is shown in [mettre article] that if the boundary of augmented accessible set is the fitted ellipse then the shooting function T_2 is the identity function. Due to the error of the approximation, the function T_2 is not the identity, but we hope that this function is close to this ideal function, and so that the convergence of T_2 is faster than the one of S_2. ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"The following code compares the convergence of these two shooting function. We also study for these functions the domain of initial guesses for the solver that hit the bounds -1 1 during the solving process. ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"## initial guesses\nN = 1000; ξ = range(-1.,1.,N) \n\n# initialization of the matrix\nfnorms_S2 = zeros(N, 100); fnorms_T2 = zeros(N, 100) \nfnorms_T2_in_x = zeros(N, 100)\niterates_S2 = zeros(N, 100); iterates_T2 = zeros(N, 100)\nconv_S2 = zeros(N,1); conv_T2 = zeros(N,1) \n#-1 : not converged; 1 converged; 0: converged but hit bounds\n\n# intermediate function to get value of S from T iterates\nT₂_(q0) = abs(q0) < 1 ? S(p₀(transpose(A)*[η.(q0),q0])) : sign(q0) * tf - xT\n\nfor i = 1:N\n # remove old iterates\n global iterate_S2 = Vector{Float64}()\n global iterate_T2 = Vector{Float64}()\n\n # solve \n q_sol_S2 = fsolve(S₂!, JS₂!, [ξ[i]], show_trace = false, tracing = true)\n q_sol_T2 = fsolve(T₂!, JT₂!, [ξ[i]], show_trace = false, tracing = true)\n\n # store results is converged \n if q_sol_S2.converged\n fnorm_S2 = [q_sol_S2.trace.trace[j].fnorm for j ∈ 1:length(q_sol_S2.trace.trace)]\n iterates_S2[i,1:length(iterate_S2)] = iterate_S2\n conv_S2[i] = length(findall(x-> abs(x) > 1, iterate_S2)) == 0\n fnorms_S2[i,1:length(fnorm_S2)] = fnorm_S2\n else\n conv_S2[i] = -1\n end\n if q_sol_T2.converged\n fnorm_T2 = [q_sol_T2.trace.trace[j].fnorm for j ∈ 1:length(q_sol_T2.trace.trace)]\n iterates_T2[i,1:length(iterate_T2)] = iterate_T2\n conv_T2[i] = length(findall(x-> abs(x) > 1, iterate_T2)) == 0\n fnorms_T2[i,1:length(fnorm_T2)] = fnorm_T2 \n fnorms_T2_in_x[i, 1:length(iterate_T2)] = abs.(T₂_.(iterate_T2))\n else\n conv_T2[i] = -1\n end\nend\n\n# mean \nmean_fnorms_S2 = mean(fnorms_S2, dims = 1)\nmean_fnorms_T2 = mean(fnorms_T2, dims = 1)\nmean_fnorms_T2_in_x = mean(fnorms_T2_in_x, dims = 1)\n\n# remove zeros with tolerance ε \nε = 1e-9;\nmean_fnorms_S2 = mean_fnorms_S2[1:findall(x -> x < ε, mean_fnorms_S2)[1][2]]\nmean_fnorms_T2 = mean_fnorms_T2[1:findall(x -> x < ε, mean_fnorms_T2)[1][2]]\nmean_fnorms_T2_in_x = mean_fnorms_T2_in_x[1:findall(x -> x < ε, mean_fnorms_T2_in_x)[1][2]]\n\n# plots\nplt_1 = plot(0:length(mean_fnorms_S2)-1, mean_fnorms_S2, label = \"S₂\")\nplot!(0:length(mean_fnorms_T2)-1, mean_fnorms_T2, label = \"T₂\")\nplot!(0:length(mean_fnorms_T2_in_x)-1, mean_fnorms_T2_in_x, label = \"T₂ in x\")\nplot!(yaxis = :log10, xlim = [0, 30], ylim = [ε, 10], xlabel = \"Error\", ylabel = \"Iterations\")\n\nplt_21 = plot(ξ, S₂, color = :black, label = \"\")\ncolor = [conv_S2[i]==1 ? :green : conv_S2[i] == 0 ? :blue : :red for i ∈ 1:N]\nscatter!(ξ, S₂, color = color, markerstrokecolor = color, marker = 2, label =\"\")\nplot!([-1,1], [xT,xT], color = :black, label = \"\")\nplot!(xlabel = \"p₀\", ylabel = \"S₂\", label = \"\")\n\nplt_22 = plot(ξ, T₂, color = :black, label = \"\")\ncolor = [conv_T2[i]==1 ? :green : conv_T2[i] == 0 ? :blue : :red for i ∈ 1:N]\nscatter!(ξ, T₂, markercolor = color, markerstrokecolor = color, marker = 2, label = \"\")\nplot!([-1,1], [yT,yT], color = :black, label = \"\")\nplot!(xlabel = \"q₀\", ylabel = \"T₂\", label = \"\")\nscatter!(1,1, markercolor = :green, markerstrokecolor = :green, marker = 2, label = \"converged\")\nscatter!(1,1, markercolor = :blue, markerstrokecolor = :blue, marker = 2, label = \"converged but hit bounds\")\nscatter!(1,1, markercolor = :red, markerstrokecolor = :red, marker = 2, label = \"not converged\")\n\nplt_2 = plot(plt_21, plt_22, layout = (1,2))\nplot(plt_1, plt_2, layout = grid(2,1, heights = [0.5, 0.5]), size=(700, 600))","category":"page"},{"location":"index.html#Optimal-control-problem","page":"Introduction","title":"Optimal control problem","text":"","category":"section"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"We consider the following optimal control problem ","category":"page"},{"location":"index.html","page":"Introduction","title":"Introduction","text":" left beginarrayll\n displaystyle min_xu int_t_0^t_f x(t) mathrm dt 1em\n textscdot x(t) = u(t) tin t_0 t_fmathrmae 05em\n phantommathrmsc u(t) in -11 tin t_0 t_f 05em\n phantommathrmsc x(t_0) = x_0 quad x(t_f) = x_f\n endarray right","category":"page"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"with x_0, t_0, x_f and t_f fixed. This problem is simple, and can be analytically solve without the use of numerical method. However, the goal is to solve this problem by indirect shooting. ","category":"page"},{"location":"index.html#Dependencies","page":"Introduction","title":"Dependencies","text":"","category":"section"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"All the numerical simulations to generate this documentation from MRI.jl are performed with the following packages.","category":"page"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"using Pkg\nPkg.status()","category":"page"},{"location":"2D-example.html#Indirect-method","page":"Indirect shooting","title":"Indirect method","text":"","category":"section"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"We introduce the pseudo-Hamiltonian ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":" h(xpp^0u) = p^0 x + p u","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"For the sake of simplicity, we consider in this notebook only the normal case, and we fix p^0 = -1. According to the Pontryagin maximum principle, the maximizing control is given by u(xp) to mathrmsign(p). This function is non-differentiable, and may lead to numerical issues. ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"using OptimalControl # hide\nusing Plots # hide\nusing ForwardDiff # hide\nusing DifferentialEquations # hide\nusing MINPACK # hide\n\nt0 = 0 # initial time\nx0 = 0 # initial state\ntf = 5 # final time\nxf = 0 # final state\n\n@def ocp begin # problem definition\n\n t ∈ [ t0, tf ], time\n x ∈ R, state\n u ∈ R, control\n\n x(t0) == x0\n x(tf) == xf\n\n ẋ(t) == u(t) \n\n ∫( x(t) ) → min\n\nend\nnothing # hide","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"using MINPACK\nfunction fsolve(f, j, x; kwargs...)\n try\n MINPACK.fsolve(f, j, x; kwargs...)\n catch e\n println(\"Erreur using MINPACK\")\n println(e)\n println(\"hybrj not supported. Replaced by hybrd even if it is not visible on the doc.\")\n MINPACK.fsolve(f, x; kwargs...)\n end\nend\nfunction fsolve(f, x; kwargs...)\n MINPACK.fsolve(f, x; kwargs...)\nend","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"Thanks to the control-toolbox, the flow varphi of the (true) Hamiltonian","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":" H(xp) = h(xp-1 u(xp)) = p^0 x + lvert p rvert ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"is given by the function textttFlow. The shooting function S colon mathbbR to mathbbR is defined by","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":" S(p_0) = pi big( varphi(t_0 x_0 p_0 t_f) big) - x_f","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"where pi (xp) = x is the classical x-space projection.","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"ϕ = Flow(ocp, (x,p) -> sign(p)) # flow with maximizing control \nπ((x,p)) = x; # projection on state space\n\nS(p0) = π( ϕ(t0, x0, p0, tf) ) - xf; # shooting function\nnle = p0 -> [S(p0[1])] # intermediate function\n\n# Plot\nplot(range(-7, 2, 500), S, xlim = [-7, 2])\nplot!([-7,2], [0,0], color = :black)\nplot!(xlabel = \"p0\", ylabel = \"S\", legend=false)","category":"page"},{"location":"2D-example.html#Finite-difference-method","page":"Indirect shooting","title":"Finite difference method","text":"","category":"section"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"The main goal now is to find the zero of S. To this purpose, we use the numerical solver texttthybrd1 given in the package textttMINPACKjl. If we don't provide the Jacobian J_S of S to the solver, the finite difference method is used to approximate it. ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"ξ = [-1.0] # initial guess\nS!(s, ξ) = (s[:] .= S(ξ[1]); nothing) # intermediate function\np0_sol = fsolve(S!, ξ, show_trace = true) # solve\nprintln(p0_sol)","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"sol = ϕ((t0, tf), x0, p0_sol.x) # get the optimal trajectory\nplot(sol) # plot","category":"page"},{"location":"2D-example.html#Automatic-differentiation-(wrong-way)","page":"Indirect shooting","title":"Automatic differentiation (wrong way)","text":"","category":"section"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"Now, we want to provide J_S to the solver, thanks to the textttForwardDiffjl package. This Jacobian is computed with the variational equation, and leads to a false result in our case.","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"Details. ( To be removed )","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"Denoting z = (xp), we have ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":" varphi(t_0 z_0 t_f) = z_0 + int_t_0^t_f vec Hbig(varphi(t_0 z_0 t)big) mathrm dt ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"If we assume that z_0 to varphi(t_0 z_0 t_f) is differentiable, we have ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":" fracpartial varphipartial z_0(t_0 z_0 t_f)cdot delta z_0 = delta z_0 + int_t_0^t_f vec Hbig(varphi(t_0 z_0 t)big)cdot left( fracpartial varphipartial z_0(t_0 z_0 t) cdot delta z_0 right) mathrm dt ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"and so, z_0 to fracpartial varphipartial z_0(t_0 z_0 t_f)cdot delta z_0 is solution of the variational equations","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":" fracpartial delta zpartial t(t) = vec Hbig(varphi(t_0 z_0 t_f)big) cdot delta z(t) qquad delta z(t_0) = delta z_0","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"In the studied optimal control problem, we have ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":" vec H(xp) = (mathrmsign(p) -1) ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"and so, we have vec H(z) = 0_2 almost everywhere, which implies","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":" fracpartial varphipartial z_0(t_0 z_0 t_f) cdot delta z_0 = mathrmexpbig((t_f-t_0) 0_2 big)cdot delta z_0 = delta z_0","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"The Jacobian of the shooting function is then given by ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":" S(p_0) = pi left( fracpartial varphipartial p_0(t_0 x_0 p_0 t_f) right) = pi left( fracpartial varphipartial z_0(t_0 z_0 t_f) cdot (01) right) = pi(01) = 0 ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"
\n
\n \n Details.\n
\n
\n

Denoting z=(x,p)z = (x,p)z=(x,p), we have

\n

φ(t0,z0,tf)=z0+∫t0tfH⃗(φ(t0,z0,t)) dt. \\varphi(t_0, z_0, t_f) = z_0 + \\int_{t_0}^{t_f} \\vec H\\big(\\varphi(t_0, z_0, t)\\big) \\,\\mathrm dt. φ(t0​,z0​,tf​)=z0​+∫t0​tf​​H(φ(t0​,z0​,t))dt.

\n

If we assume that z0→φ(t0,z0,tf)z_0 \\to \\varphi(t_0, z_0, t_f)z0​→φ(t0​,z0​,tf​) is differentiable, we have

\n

∂φ∂z0(t0,z0,tf)⋅δz0=δz0+∫t0tfH⃗′(φ(t0,z0,t))⋅(∂φ∂z0(t0,z0,t)⋅δz0) dt, \\frac{\\partial \\varphi}{\\partial z_0}(t_0, z_0, t_f)\\cdot \\delta z_0 = \\delta z_0 + \\int_{t_0}^{t_f} \\vec H'\\big(\\varphi(t_0, z_0, t)\\big)\\cdot \\left( \\frac{\\partial \\varphi}{\\partial z_0}(t_0, z_0, t) \\cdot \\delta z_0 \\right) \\,\\mathrm dt, ∂z0​∂φ​(t0​,z0​,tf​)⋅δz0​=δz0​+∫t0​tf​​H′(φ(t0​,z0​,t))⋅(∂z0​∂φ​(t0​,z0​,t)⋅δz0​)dt,

\n

and so, z0→∂φ∂z0(t0,z0,tf)⋅δz0z_0 \\to \\frac{\\partial \\varphi}{\\partial z_0}(t_0, z_0, t_f)\\cdot \\delta z_0z0​→∂z0​∂φ​(t0​,z0​,tf​)⋅δz0​ is solution of the variational equations

\n

∂δz∂t(t)=H⃗′(φ(t0,z0,tf))⋅δz(t),δz(t0)=δz0. \\frac{\\partial \\delta z}{\\partial t}(t) = \\vec H'\\big(\\varphi(t_0, z_0, t_f)\\big) \\cdot \\delta z(t), \\qquad \\delta z(t_0) = \\delta z_0.∂t∂δz​(t)=H′(φ(t0​,z0​,tf​))⋅δz(t),δz(t0​)=δz0​.

\n

In the studied optimal control problem, we have

\n

H⃗(x,p)=(sign(p),−1) \\vec H(x,p) = (\\mathrm{sign}(p), -1) H(x,p)=(sign(p),−1)

\n

and so, we have H⃗′(z)=02\\vec H'(z) = 0_2H′(z)=02​ almost everywhere, which implies

\n

∂φ∂z0(t0,z0,tf)⋅δz0=exp((tf−t0)02)⋅δz0=δz0. \\frac{\\partial \\varphi}{\\partial z_0}(t_0, z_0, t_f) \\cdot \\delta z_0 = \\mathrm{exp}\\big((t_f-t_0) 0_2 \\big)\\cdot \\delta z_0 = \\delta z_0.∂z0​∂φ​(t0​,z0​,tf​)⋅δz0​=exp((tf​−t0​)02​)⋅δz0​=δz0​.

\n

The Jacobian of the shooting function is then given by

\n

S′(p0)=π(∂φ∂p0(t0,x0,p0,tf))=π(∂φ∂z0(t0,z0,tf)⋅(0,1))=π(0,1)=0. S'(p_0) = \\pi \\left( \\frac{\\partial \\varphi}{\\partial p_0}(t_0, x_0, p_0, t_f) \\right) = \\pi \\left( \\frac{\\partial \\varphi}{\\partial z_0}(t_0, z_0, t_f) \\cdot (0,1) \\right) = \\pi(0,1) = 0. S′(p0​)=π(∂p0​∂φ​(t0​,x0​,p0​,tf​))=π(∂z0​∂φ​(t0​,z0​,tf​)⋅(0,1))=π(0,1)=0.

\n
\n
","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"ξ = [-1.0] # initial guess\nJS(ξ) = ForwardDiff.jacobian(p -> [S(p[1])], ξ) # compute jacobian by forward differentiation\nprintln(\"ξ = \", ξ[1])\nprintln(\"JS(ξ) : \", JS(ξ)[1])","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"However, the solver texttthybrd1 uses rank 1 approximation to actualize the Jacobian insted of compute it at each iteration, which imply that it still converges to the solution even if the given Jacobian is completely false.","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"JS!(js, ξ) = (js[:] .= JS(ξ); nothing) # intermediate function\np0_sol = fsolve(S!, JS!, ξ, show_trace = true) # solve\nprintln(p0_sol)","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"sol = ϕ((t0, tf), x0, p0_sol.x) # get the optimal trajectory\nplt = plot(sol) # plot","category":"page"},{"location":"2D-example.html#Automatic-differentiation-(good-way)","page":"Indirect shooting","title":"Automatic differentiation (good way)","text":"","category":"section"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"The goal is to provide the true Jacobian of S by using the textttForwardDiff package, and so we need to indicate to the solver that the dynamic of the system change when p = 0.","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"To understang why we need to give this information to the solver, see the following details. ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"Details. (To be removed)","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"The problem is that the Hamiltonian H is not differentiable everywhere due to the maximizing control. This control is bang-bang (u = 1 and u = -1). ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"Let now construct the two smooth Hamiltonians associated to these two controls ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":" H^+(xp) = h(xp-11) = -x + p qquad textand qquad H^-(xp) = h(xp-1-1) = -x - p ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"Their associated vector fields are given by ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":" vec H^+(xp) = (11) qquad textand qquad vec H^-(xp) = (-1 1) ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"and their associated flow correspond to ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":" varphi^+(t_0 z_0 t_f) = z_0 + left( beginarrayc 1 1 endarray right) (t_f -t_0)\n qquad textand qquad \n varphi^-(t_0 z_0 t_f) = z_0 + left( beginarrayc -1 phantom- 1 endarray right) (t_f -t_0)","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"If we assume that the optimal structure of the problem is negative then positive bangs, then the associated flow is defined by ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":" varphi(t_0 z_0 t_f) = varphi^+ big( t_1(z_0) varphi^-big(t_0 z_0 t_1(z_0)big) t_f big)","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"with the following condition ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":" pi_p big( varphi^-(t_0 z_0 t_1(z_0)) big) = 0","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"where pi_p(xp) = p is the classical p-space projection. By devlopping this last condition, an explicit form of the function t_1(cdot) is given by ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":" t_1(x_0 p_0) = t_0 - p_0","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"Finally, we have ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":" beginalign*\n fracpartial varphipartial z_0 \n = fracpartial varphi^+partial t_0 fracpartial t_1partial z_0 + fracvarphi^+partial z_0 left( fracpartial varphi^-partial z_0 + fracpartial varphi^-partial t_f fracpartial t_1partial z_0 right) \n = left( beginarrayc -1 -1 endarray right) left( beginarraycc0 -1 endarray right) + left( beginarraycc 1 0 0 1 endarray right) left\n left( beginarraycc 1 0 0 1 endarray right) + left( beginarrayc -1 phantom - 1 endarray right) left( beginarraycc0 -1 endarray right)\n right \n = left( beginarraycc 0 1 0 1 endarray right) + left( beginarraycc 1 0 0 1 endarray right) + left( beginarraycc 0 phantom -1 0 -1 endarray right) \n = left( beginarraycc 1 2 0 1 endarray right)\n endalign*","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"and so, we have that ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":" S(p_0) = pi left( fracpartial varphipartial p_0(t_0 x_0 p_0 t_f) right) = pi left( fracpartial varphipartial z_0(t_0 z_0 t_f) cdot (01) right) = pi(21) = 2","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"
\n
\n \n Details.\n
\n
\n

The problem is that the Hamiltonian HHH is not differentiable everywhere due to the maximizing control. This control is bang-bang (u=1u = 1u=1 and u=−1u = -1u=−1).

\n

Let now construct the two smooth Hamiltonians associated to these two controls

\n

H+(x,p)=h(x,p,−1,1)=−x+pandH−(x,p)=h(x,p,−1,−1)=−x−p. H^+(x,p) = h(x,p,-1,1) = -x + p \\qquad \\text{and} \\qquad H^-(x,p) = h(x,p,-1,-1) = -x - p. H+(x,p)=h(x,p,−1,1)=−x+pandH−(x,p)=h(x,p,−1,−1)=−x−p.

\n

Their associated vector fields are given by

\n

H⃗+(x,p)=(1,1)andH⃗−(x,p)=(−1,1), \\vec H^+(x,p) = (1,1) \\qquad \\text{and} \\qquad \\vec H^-(x,p) = (-1, 1), H+(x,p)=(1,1)andH−(x,p)=(−1,1),

\n

and their associated flow correspond to

\n

φ+(t0,z0,tf)=z0+(11)(tf−t0)andφ−(t0,z0,tf)=z0+(−1−1)(tf−t0). \\varphi^+(t_0, z_0, t_f) = z_0 + \\left( \\begin{array}{c} 1 \\\\ 1 \\end{array} \\right) (t_f -t_0)\n \\qquad \\text{and} \\qquad \n \\varphi^-(t_0, z_0, t_f) = z_0 + \\left( \\begin{array}{c} -1 \\\\ \\phantom{-} 1 \\end{array} \\right) (t_f -t_0).φ+(t0​,z0​,tf​)=z0​+(11​)(tf​−t0​)andφ−(t0​,z0​,tf​)=z0​+(−1−1​)(tf​−t0​).

\n

If we assume that the optimal structure of the problem is negative then positive bangs, then the associated flow is defined by

\n

φ(t0,z0,tf)=φ+(t1(z0),φ−(t0,z0,t1(z0)),tf), \\varphi(t_0, z_0, t_f) = \\varphi^+ \\big( t_1(z_0), \\varphi^-\\big(t_0, z_0, t_1(z_0)\\big), t_f \\big),φ(t0​,z0​,tf​)=φ+(t1​(z0​),φ−(t0​,z0​,t1​(z0​)),tf​),

\n

with the following condition

\n

πp(φ−(t0,z0,t1(z0)))=0, \\pi_p \\big( \\varphi^-(t_0, z_0, t_1(z_0)) \\big) = 0,πp​(φ−(t0​,z0​,t1​(z0​)))=0,

\n

where πp(x,p)=p\\pi_p(x,p) = pπp​(x,p)=p is the classical ppp-space projection. By devlopping this last condition, an explicit form of the function t1(⋅)t_1(\\cdot)t1​(⋅) is given by

\n

t1(x0,p0)=t0−p0. t_1(x_0, p_0) = t_0 - p_0.t1​(x0​,p0​)=t0​−p0​.

\n

Finally, we have

\n

∂φ∂z0=∂φ+∂t0∂t1∂z0+φ+∂z0(∂φ−∂z0+∂φ−∂tf∂t1∂z0)=(−1−1)(0−1)+(1001)[(1001)+(−1−1)(0−1)]=(0101)+(1001)+(0−10−1)=(1201) \\begin{align*}\n \\frac{\\partial \\varphi}{\\partial z_0} \n &= \\frac{\\partial \\varphi^+}{\\partial t_0} \\frac{\\partial t_1}{\\partial z_0} + \\frac{\\varphi^+}{\\partial z_0} \\left( \\frac{\\partial \\varphi^-}{\\partial z_0} + \\frac{\\partial \\varphi^-}{\\partial t_f} \\frac{\\partial t_1}{\\partial z_0} \\right) \\\\\n &= \\left( \\begin{array}{c} -1 \\\\ -1 \\end{array} \\right) \\left( \\begin{array}{cc}0 & -1 \\end{array} \\right) + \\left( \\begin{array}{cc} 1 & 0 \\\\ 0 & 1 \\end{array} \\right) \\left[\n \\left( \\begin{array}{cc} 1 & 0 \\\\ 0 & 1 \\end{array} \\right) + \\left( \\begin{array}{c} -1 \\\\ \\phantom - 1 \\end{array} \\right) \\left( \\begin{array}{cc}0 & -1 \\end{array} \\right)\n \\right] \\\\\n &= \\left( \\begin{array}{cc} 0 & 1 \\\\ 0 & 1 \\end{array} \\right) + \\left( \\begin{array}{cc} 1 & 0 \\\\ 0 & 1 \\end{array} \\right) + \\left( \\begin{array}{cc} 0 & \\phantom -1 \\\\ 0 & -1 \\end{array} \\right) \\\\ \n &= \\left( \\begin{array}{cc} 1 & 2 \\\\ 0 & 1 \\end{array} \\right)\n \\end{align*}∂z0​∂φ​​=∂t0​∂φ+​∂z0​∂t1​​+∂z0​φ+​(∂z0​∂φ−​+∂tf​∂φ−​∂z0​∂t1​​)=(−1−1​)(0​−1​)+(10​01​)[(10​01​)+(−1−1​)(0​−1​)]=(00​11​)+(10​01​)+(00​−1−1​)=(10​21​)​

\n

and so, we have that

\n

S′(p0)=π(∂φ∂p0(t0,x0,p0,tf))=π(∂φ∂z0(t0,z0,tf)⋅(0,1))=π(2,1)=2. S'(p_0) = \\pi \\left( \\frac{\\partial \\varphi}{\\partial p_0}(t_0, x_0, p_0, t_f) \\right) = \\pi \\left( \\frac{\\partial \\varphi}{\\partial z_0}(t_0, z_0, t_f) \\cdot (0,1) \\right) = \\pi(2,1) = 2.S′(p0​)=π(∂p0​∂φ​(t0​,x0​,p0​,tf​))=π(∂z0​∂φ​(t0​,z0​,tf​)⋅(0,1))=π(2,1)=2.

\n
\n
","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"To provide this change of dynamic to the solver, we need to use a callback during the integration that will execute the function textttaffect when textttcondition(xp) = 0.","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"For us, the condition is given by (xp) to p. For the textttaffect function, we use a global parameter alpha. This parameter will be set to pm 1 at the beginning of the integration and it sign will change with the textttaffect function. ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"Thanks to the textttcontrol-toolbox package, the created callback can be easily pass to the integrator throught the textttFlow function.","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"global α # parameter: ̇p(t) = α with α = ±1\n\nfunction condition(z, t, integrator) # event when condition(x,p) == 0\n x,p = z\n return p\nend\n\nfunction affect!(integrator) # action when condition == 0 \n global α = -α\n nothing\nend\n\ncb = ContinuousCallback(condition, affect!) # callback \n\nφ_ = Flow(ocp, (x,p) -> α, callback = cb) # intermediate flow\n\nfunction φ(t0, x0, p0, tf; kwargs...) # flow\n global α = sign(p0)\n return φ_(t0, x0, p0, tf; kwargs...)\nend\n\nfunction φ((t0, tf), x0, p0; kwargs...) # flow for plot\n global α = sign(p0)\n return φ_((t0, tf), x0, p0; kwargs...)\nend\n\nShoot(p0) = π( φ(t0, x0, p0, tf) ) - xf # shooting function\n\nξ = [-1.0] # initial guess\nJShoot(ξ) = ForwardDiff.jacobian(p -> [Shoot(p[1])], ξ) # compute jacobian by forward differentiation\nprintln(\"ξ = \", ξ[1])\nprintln(\"JS(ξ) : \", JShoot(ξ)[1])\nShoot!(shoot, ξ) = (shoot[:] .= Shoot(ξ[1]); nothing) # intermediate function\nJShoot!(jshoot, ξ) = (jshoot[:] .= JShoot(ξ); nothing) # intermediate function\n\np0_sol = fsolve(Shoot!, JShoot!, ξ, show_trace = true) # solve\nprintln(p0_sol)","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"# get optimal trajectory\nsol_ = φ((t0, tf), x0, p0_sol.x[1], saveat=range(t0, tf, 500)) \n\n# plot\nsol = OptimalControl.OptimalControlSolution(sol_)\nt = sol.times\nx = sol.state\np = sol.costate\nu = sign ∘ p\n\nplt_x = plot(t, x, label = \"x\")\nplt_p = plot(t, p, label = \"p\")\nplt_u = plot(t, u, label = \"u\")\n\nplt_xp = plot(plt_x, plt_p, layout=(1, 2))\nplot(plt_xp, plt_u, layout = (2, 1))","category":"page"}] +[{"location":"2D-preconditioner.html#Augmented-formulation","page":"Geometric preconditioner","title":"Augmented formulation","text":"","category":"section"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"We consider the augmented formulation of the same problem as before","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":" left beginarrayll\n displaystyle min_xu x^0(t_f) 1em\n textscdot x^0(t) = x_1(t) tin t_0 t_fmathrmae 05em\n phantommathrmsc dot x_1(t) = u(t) tin t_0 t_fmathrmae 05em\n phantommathrmsc u(t) in -11 tin t_0 t_f 05em\n phantommathrmsc x(t_0) = x_0 quad x_1(t_f) = x_f\n endarray right","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"with fixed x_0, t_0, x_f and t_f, and where x = (x^0 x_1). The augmented Hamiltonian of this problem is given by ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":" H(xp) = p^0 x_1 + lvert p_1 rvert","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"where p = (p^0 p_1) is the augmented costate, composed by the costate p_1 associated to the state x_1 and multiplier p^0 associated to the cost x^0 which is no longer fixed to -1. The flow varphi associated to the true Hamiltonian H is computed thanks to the textttflow and the textttHamiltonian functions.","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"The general shooting function S colon mathbb R^2 to mathbb R associated to this Hamiltonian is defined by ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":" S(p_0) = pi_x big( varphi(t_0 x_0 p_0 t_f) big)","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"where pi_x(x^0 x_1 p^0 x_1) = x_1 is still the state projection. ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"We are interested in two normalization of this shooting function, thanks to the homogeneity of the BC-extremals on the augmented costate, which are ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":" S_1(p_0) = S(-1 p_0) qquad textand qquad S_2(p_0) = S big( eta(p) p_0 big)","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"where the function eta colon -1 1 to mathbb R is defined by ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":" eta(p) = -sqrt 1 - p^2","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"The following code creates and plots these shooting functions. ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"using OptimalControl\nusing Plots\nusing ForwardDiff\nusing DifferentialEquations\nusing MINPACK\nusing Statistics\n\nglobal α = 1\n\nfunction condition(z, t, integrator) # event when condition(z,t,integrator) == 0\n x⁰,x,p⁰,p = z\n return p\nend\n\nfunction affect!(integrator) # action when condition == 0 \n global α = -α\n nothing\nend\n\ncb = ContinuousCallback(condition, affect!) # callback \nH(x, p) = p[1] * x[2] + α * p[2]\nϕ_ = Flow(Hamiltonian(H), callback = cb) # flow with maximizing control \n\nfunction ϕ(t0, x0, p0, tf; kwargs...) \n if p0[2] == 0\n global α = -sign(p0[1])\n else\n global α = sign(p0[2])\n end\n return ϕ_(t0, x0, p0, tf; kwargs...)\nend\n\nfunction ϕ((t0, tf), x0, p0; kwargs...) # flow for plot\n if p0[2] == 0\n global α = -sign(p0[1])\n else\n global α = sign(p0[2])\n end\n return ϕ_((t0, tf), x0, p0; kwargs...)\nend\n\nt0 = 0 # initial time\nx0 = [0,0] # initial augmented state\ntf = 5 # final time\nxT = 0 # final state\n\nπ((x,p)) = x[2] # projection on state space\nη(p0) = -sqrt.(1 - p0.^2) # function η(⋅)\n\nS(p0) = π( ϕ(t0, x0, p0, tf) ) - xT # general hooting function\n\nS₁(p0) = S([-1, p0]) # normalization 1\nS₂(p0) = abs(p0) < 1 ? S([η(p0), p0]) : sign(p0)*tf - xT # normalization 2\n\n# Plot\nplt_S1 = plot(range(-7, 2, 500), S₁ , color = 2, label = \"\")\nplot!(xlabel = \"p₀\", ylabel = \"S₁\", xlim = [-7,2])\n\nplt_S2 = plot(range(-1, 1, 500), S₂, color = 3)\nplot!(xlabel = \"p₀\", ylabel = \"S₂\", xlim = [-1,1], legend=false)\n\nS_(p⁰, p) = S([p⁰, p])\nplt_S = surface(range(0, -1.5, 100), range(-8, 2, 100), S_, camera = (30,30))\nsurface!(xlabel = \"p⁰\", ylabel = \"p₀\", zlabel = \"φₓ\", xflip = true)\nplot3d!(-1*ones(100), range(-8, 2, 100), S₁.(range(-8, 2, 100)), label = \"S₁\")\nplot3d!(η.(range(-1, 1, 100)), range(-1, 1, 100), S₂.(range(-1, 1, 100)), label = \"S₂\")\n\nplt_S12 = plot(plt_S1, plt_S2, layout = (1,2))\nplt_total = plot(plt_S, plt_S12, layout = grid(2,1, heights = [2/3, 1/3]), size=(700, 600))\n\n@gif for i ∈ [range(30, 90, 50); 90*ones(25); range(90, 30, 50); 30*ones(25)]\n plot!(plt_total[1], camera=(i,i), \n zticks = i==90 ? false : true,\n zlabel = i==90 ? \"\" : \"S\" )\nend","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"As done before, we use the solver texttthybrd1 from the textttMINPACKjl package to find a zero of S_2.","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"global iterate_S2 = Vector{Float64}() # global vector to store iterates of the solver\n\nfunction S₂!(s₂, ξ) # intermediate function\n push!(iterate_S2, ξ[1]) # stock the iterate\n return (s₂[:] .= S₂(ξ[1]); nothing) \nend\n\nJS₂(ξ) = ForwardDiff.jacobian(p0 -> [S₂(p0[1])], ξ) # compute jacobian by forward differentiation\nJS₂!(js₂, ξ) = (js₂[:] .= JS₂(ξ); nothing) # intermediate function\n\nξ = [-0.5] # initial guess\np0_sol = fsolve(S₂!, JS₂!, ξ, show_trace = true) # solve\nprintln(p0_sol) # print solution","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"sol = ϕ((t0, tf), x0, [η(p0_sol.x[1]), p0_sol.x[1]], \n saveat=range(t0, tf, 500)) # get optimal trajectory\n\n# plot\nt = sol.t\nx⁰ = [sol.u[i][1] for i in 1:length(sol.u)]\nx = [sol.u[i][2] for i in 1:length(sol.u)]\np⁰ = [sol.u[i][3] for i in 1:length(sol.u)]\np = [sol.u[i][4] for i in 1:length(sol.u)]\nu = sign.(p)\n\nplt_x⁰ = plot(t, x⁰, label = \"x⁰\")\nplt_x = plot(t, x , label = \"x\" )\nplt_p⁰ = plot(t, p⁰, label = \"p⁰\", ylim=[-1,1])\nplt_p = plot(t, p , label = \"p\" )\nplt_u = plot(t, u, label = \"u\" )\n\nplt_xp = plot(plt_x⁰, plt_p⁰, plt_x, plt_p, layout=(2, 2))\nplot(plt_xp, plt_u, layout = grid(2,1, heights = [2/3, 1/3]), size=(700, 600))","category":"page"},{"location":"2D-preconditioner.html#Construction-of-the-geometric-preconditioner","page":"Geometric preconditioner","title":"Construction of the geometric preconditioner","text":"","category":"section"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"The goal is now to use the geometric preconditioning method proposed in [mettre article]. For this purpose, the first step is to create some points in the boundary of the accessible augmented set, and to fit an ellipse on these points. ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"The second step is to create the linear diffeomorphism phi colon mathbb R^2 to mathbb R^2 (x^0x_1) to Ax + B that transforms the fitted ellipse into the unit circle, and that satisfy the condition ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":" fracpartial phipartial x^0 = k e_1 ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"with k0 and e_1 = (10). In this context, we denote ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":" A = left( beginarraycc k A_x^0 0 A_x_1 endarray right)\n qquad textand qquad \n B = left( beginarrayc B_x^0 B_x_1 endarray right)","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"This diffeomorphism is given from the semi-axis ab 0, the angle theta in 0 fracpi2 between the semi-axis b and the x-axis, and the center c in mathbb R^2 by ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":" phi(x) = r(-beta_0) s(a^-1 b^-1) r(theta) (x - c)","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"where r and s correspond respectively to the rotation and the scale matrix, defined by ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":" r(theta) = left( beginarraycc phantom - cos(theta) sin(theta) -sin(theta) cos(theta) endarray right)\n qquad textand qquad\n s(ab) = left( beginarraycc a 0 0 b endarray right)","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"and where beta_0 = arctan left(fraca sin(theta)b cos(theta) right). ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"
\n
\n \n fit_ellipse — Function\n
\n
function fit_ellipse(x, y)\n    M = hcat(x.^2, x.*y, y.^2, x, y)                            # quadratic form of ellipse\n    p = M\\ones(length(x))                                       # fit parameters for the ellipse\n    A, B, C, D, E = p\n    F = -1.0\n    # calculate the parameters from quadratic\n    Δ = B^2 - 4*A*C\n    Λ = (A-C)^2 + B^2\n    b, a = [-sqrt(clamp( 2*(A*E^2 + C*D^2 - B*D*E + Δ*F)*\n            ( (A+C) + op(sqrt(Λ)) ), 0, Inf)) / Δ   for op in (+, -)]\n    θ = atan(-B, C-A)/2\n    c = [(2*C*D - B*E)/Δ, (2*A*E - B*D)/Δ]\n    return a, b, -θ+Base.π/2, c\nend
\n
\n
","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"function fit_ellipse(x, y) # hide\n M = hcat(x.^2, x.*y, y.^2, x, y) # hide\n p = M\\ones(length(x)) # hide\n A, B, C, D, E = p # hide\n F = -1.0 # hide\n Δ = B^2 - 4*A*C # hide\n Λ = (A-C)^2 + B^2 # hide\n b, a = [-sqrt(clamp( 2*(A*E^2 + C*D^2 - B*D*E + Δ*F) * ( (A+C) + op(sqrt(Λ)) ), 0, Inf)) / Δ for op in (+, -)] # hide\n θ = atan(-B, C-A)/2 # hide\n c = [(2*C*D - B*E)/Δ, (2*A*E - B*D)/Δ] # hide\n return a, b, -θ+Base.π/2, c # hide\nend # hide\nnothing # hide","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"n = 15 # number of points for fit : 2n\nn_ = 100 # number of points for plot: 2n_\np0 = [[[-1, i] for i ∈ range(-tf, 0, n)]; # initial costate for fit\n [[1, i] for i ∈ range(tf, 0, n)]] \np0_ = [[[-1, i] for i ∈ range(-tf, 0, n_)]; # initial costate for plot\n [[1, i] for i ∈ range(tf, 0, n_)]]\n\nx = zeros(2, 2*n); p = zeros(2,2*n) # init final state and costate\nx_ = zeros(2, 2*n_); p_ = zeros(2, 2*n_)\nfor i = 1:length(p0)\n x[:,i], p[:,i] = ϕ(t0, x0, p0[i], tf) # compute flow for fit\nend\nfor i = 1:length(p0_)\n x_[:,i], p_[:,i] = ϕ(t0, x0, p0_[i], tf) # compute flow for plot\nend\n\na, b, θ, c = fit_ellipse(x[1,:], x[2,:]) # fit ellipse\nr(β) = [[cos(β), sin(β)] [-sin(β), cos(β)]] # 2x2 rotation matrix\ns(a,b) = [[a,0] [0,b]] # 2x2 scale matrix\nβ = range(-Base.π, Base.π; length = 100) # angle for plot ellipse\nxₑ = r(-θ)*s(a,b)*\n transpose(reduce(hcat,[sin.(β), cos.(β)])).+c # points of the ellipse\n\n# construction of the linear diffeomorphism φ(x) = Ax + B\nd = (a*sin(θ))/(b*cos(θ)); β₀ = atan(d) # intermediate values\nA = r(-β₀)*s(1/a,1/b)*r(θ); B = -A*c # calculate A and B\nφ(x) = A*x .+ B # function φ\n\ny = φ(x); y_ = φ(x_); yₑ = φ(xₑ) # compute φ on x, x_ and xₑ\n\n# plot\nplt_x = plot(x_[2,:], x_[1,:], label = \"\")\nscatter!(x[2,:], x[1,:], label=\"Observations\", legend = :topleft)\nplot!(xₑ[2,:], xₑ[1,:], label = \"Fitted ellipse\")\nplot!(xlim = [-15,15], ylim = [-15,15], xlabel = \"x\", ylabel = \"x⁰\")\n\nplt_y = plot(y_[2,:], y_[1,:], label = \"\")\nscatter!(y[2,:], y[1,:], label=\"\")\nplot!(yₑ[2,:], yₑ[1,:], label = \"\")\nplot!( xlabel = \"y\", ylabel = \"y⁰\")\n\nplot(plt_x, plt_y, layout = (1,2), size=(800, 400))","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"The general shooting function T colon mathbb R^2 to mathbb R in the new coordinate is defined by ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":" T(q) = A_x_1 varphi big(t_0 x_0 p_0(A^top q) t_f) + B_x_1 - y_T","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"where the function p_0 colon mathbb R^2 to mathbb R^2 correspond to the mapping between the final and the initial augmented costate, and is given by","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":" p_0(p^0 p_1) = (p^0 p_1 + p^0 t_f)","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"and y_T = A_x_1 x_T + B_x to the target in the new system of coordinates. By using the definition of S and y_T, we obtain ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"beginalign*\n T(q) = A_x_1 varphi big(t_0 x_0 p_0(A^top q) t_f) + B_x_1 - (A_x_1 x_T + B_x_1) \n = A_x_1 (S circ p_0)(A^top q)\nendalign*","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"which highlight that the proposed geometric preconditioning method is a left and right side preconditioner of the shooting function. Finally, we define the two shooting functions T_1 and T_2 by using the two method of normalization used before for the function S ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":" T_1(p) = T(-1 p) qquad textand qquad T_2(p) = T big(eta(p) p big)","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"p₀(p) = [p[1], p[2] + p[1]*tf] # function p₀(⋅)\nAₓ = A[2,2]; Bₓ = B[2] # Aₓ and Bₓ\nyT = Aₓ*xT + Bₓ # target in the new coordinate \nT(q) = Aₓ*(S ∘ p₀)(transpose(A)*q) # state flow in the new coordinates\nT₁(q) = T([-1, q]) # normalization 1\nT₂(q) = abs(q) < 1 ? T([η(q), q]) : sign(q)*(Aₓ*tf + Bₓ) # normalization 2\n\n# Plot\nplt_S1 = plot(range(-3, 3, 500), T₁ , color = 2)\nplot!([-3,3], [yT,yT], color = :black)\nplot!(xlabel = \"q₀\", ylabel = \"T₁\", xlim = [-3,3], legend=false)\n\nplt_S2 = plot(range(-1, 1, 500), T₂, color = 3)\nplot!([-1,1], [yT,yT], color = :black)\nplot!(xlabel = \"q₀\", ylabel = \"T₂\", xlim = [-1,1], legend=false)\n\nT_(q⁰, q) = T([q⁰, q])\nplt_S = surface(range(0, -1.5, 100), range(-3, 3, 100), T_, camera = (30,30))\nsurface!(xlabel = \"q⁰\", ylabel = \"q₀\", zlabel = \"T\", xflip = true)\nplot3d!(-1*ones(100), range(-3, 3, 100), T₁.(range(-3, 3, 100)), label = \"T₁\")\nplot3d!(η.(range(-1, 1, 100)), range(-1, 1, 100), T₂.(range(-1, 1, 100)), label = \"T₂\")\n\nplt_S12 = plot(plt_S1, plt_S2, layout = (1,2))\nplt_total = plot(plt_S, plt_S12, layout = grid(2,1, heights = [2/3, 1/3]), size=(700, 600))\n\n@gif for i ∈ [range(30, 90, 50); 90*ones(25); range(90, 30, 50); 30*ones(25)]\n plot!(plt_total[1], camera=(i,i), \n zticks = i==90 ? false : true,\n zlabel = i==90 ? \"\" : \"T\" )\nend","category":"page"},{"location":"2D-preconditioner.html#Preconditioned-shoot","page":"Geometric preconditioner","title":"Preconditioned shoot","text":"","category":"section"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"global iterate_T2 = Vector{Float64}() # global vector to store iterates of the solver\n\nfunction T₂!(t₂, ξ) # intermediate function\n push!(iterate_T2, ξ[1]) \n return (t₂[:] .= T₂(ξ[1]); nothing) \nend\n\nJT₂(ξ) = ForwardDiff.jacobian(q0 -> [T₂(q0[1])], ξ) # compute jacobian by forward differentiation\nJT₂!(jt₂, ξ) = (jt₂[:] .= JT₂(ξ); nothing) # intermediate function\n\nξ = [0.5] # initial guess\nq_sol = fsolve(T₂!, JT₂!, ξ, show_trace = true) # solve\nprintln(q_sol) # print solution","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"p0_sol = p₀(transpose(A)*[η(q_sol.x[1]), q_sol.x[1]]) # get the optimal initial costate in old coordinates\nsol = ϕ((t0, tf), x0, p0_sol, saveat=range(t0, tf, 500)) # get optimal trajectory\n\n# plot\nt = sol.t\nx⁰ = [sol.u[i][1] for i in 1:length(sol.u)]\nx = [sol.u[i][2] for i in 1:length(sol.u)]\np⁰ = [sol.u[i][3] for i in 1:length(sol.u)]\np = [sol.u[i][4] for i in 1:length(sol.u)]\nu = sign.(p)\n\nplt_x⁰ = plot(t, x⁰, label = \"x⁰\")\nplt_x = plot(t, x , label = \"x\" )\nplt_p⁰ = plot(t, p⁰, label = \"p⁰\", ylim=[-1,1])\nplt_p = plot(t, p , label = \"p\" )\nplt_u = plot(t, u, label = \"u\" )\n\nplt_xp = plot(plt_x⁰, plt_p⁰, plt_x, plt_p, layout=(2, 2))\nplot(plt_xp, plt_u, layout = grid(2,1, heights = [2/3, 1/3]), size=(700, 600))","category":"page"},{"location":"2D-preconditioner.html#Comparison","page":"Geometric preconditioner","title":"Comparison","text":"","category":"section"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"It is shown in [mettre article] that if the boundary of augmented accessible set is the fitted ellipse then the shooting function T_2 is the identity function. Due to the error of the approximation, the function T_2 is not the identity, but we hope that this function is close to this ideal function, and so that the convergence of T_2 is faster than the one of S_2. ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"The following code compares the convergence of these two shooting function. We also study for these functions the domain of initial guesses for the solver that hit the bounds -1 1 during the solving process. ","category":"page"},{"location":"2D-preconditioner.html","page":"Geometric preconditioner","title":"Geometric preconditioner","text":"## initial guesses\nN = 1000; ξ = range(-1.,1.,N) \n\n# initialization of the matrix\nfnorms_S2 = zeros(N, 100); fnorms_T2 = zeros(N, 100) \nfnorms_T2_in_x = zeros(N, 100)\niterates_S2 = zeros(N, 100); iterates_T2 = zeros(N, 100)\nconv_S2 = zeros(N,1); conv_T2 = zeros(N,1) \n#-1 : not converged; 1 converged; 0: converged but hit bounds\n\n# intermediate function to get value of S from T iterates\nT₂_(q0) = abs(q0) < 1 ? S(p₀(transpose(A)*[η.(q0),q0])) : sign(q0) * tf - xT\n\nfor i = 1:N\n # remove old iterates\n global iterate_S2 = Vector{Float64}()\n global iterate_T2 = Vector{Float64}()\n\n # solve \n q_sol_S2 = fsolve(S₂!, JS₂!, [ξ[i]], show_trace = false, tracing = true)\n q_sol_T2 = fsolve(T₂!, JT₂!, [ξ[i]], show_trace = false, tracing = true)\n\n # store results is converged \n if q_sol_S2.converged\n fnorm_S2 = [q_sol_S2.trace.trace[j].fnorm for j ∈ 1:length(q_sol_S2.trace.trace)]\n iterates_S2[i,1:length(iterate_S2)] = iterate_S2\n conv_S2[i] = length(findall(x-> abs(x) > 1, iterate_S2)) == 0\n fnorms_S2[i,1:length(fnorm_S2)] = fnorm_S2\n else\n conv_S2[i] = -1\n end\n if q_sol_T2.converged\n fnorm_T2 = [q_sol_T2.trace.trace[j].fnorm for j ∈ 1:length(q_sol_T2.trace.trace)]\n iterates_T2[i,1:length(iterate_T2)] = iterate_T2\n conv_T2[i] = length(findall(x-> abs(x) > 1, iterate_T2)) == 0\n fnorms_T2[i,1:length(fnorm_T2)] = fnorm_T2 \n fnorms_T2_in_x[i, 1:length(iterate_T2)] = abs.(T₂_.(iterate_T2))\n else\n conv_T2[i] = -1\n end\nend\n\n# mean \nmean_fnorms_S2 = mean(fnorms_S2, dims = 1)\nmean_fnorms_T2 = mean(fnorms_T2, dims = 1)\nmean_fnorms_T2_in_x = mean(fnorms_T2_in_x, dims = 1)\n\n# remove zeros with tolerance ε \nε = 1e-9;\nmean_fnorms_S2 = mean_fnorms_S2[1:findall(x -> x < ε, mean_fnorms_S2)[1][2]]\nmean_fnorms_T2 = mean_fnorms_T2[1:findall(x -> x < ε, mean_fnorms_T2)[1][2]]\nmean_fnorms_T2_in_x = mean_fnorms_T2_in_x[1:findall(x -> x < ε, mean_fnorms_T2_in_x)[1][2]]\n\n# plots\nplt_1 = plot(0:length(mean_fnorms_S2)-1, mean_fnorms_S2, label = \"S₂\")\nplot!(0:length(mean_fnorms_T2)-1, mean_fnorms_T2, label = \"T₂\")\nplot!(0:length(mean_fnorms_T2_in_x)-1, mean_fnorms_T2_in_x, label = \"T₂ in x\")\nplot!(yaxis = :log10, xlim = [0, 30], ylim = [ε, 10], xlabel = \"Error\", ylabel = \"Iterations\")\n\nplt_21 = plot(ξ, S₂, color = :black, label = \"\")\ncolor = [conv_S2[i]==1 ? :green : conv_S2[i] == 0 ? :blue : :red for i ∈ 1:N]\nscatter!(ξ, S₂, color = color, markerstrokecolor = color, marker = 2, label =\"\")\nplot!([-1,1], [xT,xT], color = :black, label = \"\")\nplot!(xlabel = \"p₀\", ylabel = \"S₂\", label = \"\")\n\nplt_22 = plot(ξ, T₂, color = :black, label = \"\")\ncolor = [conv_T2[i]==1 ? :green : conv_T2[i] == 0 ? :blue : :red for i ∈ 1:N]\nscatter!(ξ, T₂, markercolor = color, markerstrokecolor = color, marker = 2, label = \"\")\nplot!([-1,1], [yT,yT], color = :black, label = \"\")\nplot!(xlabel = \"q₀\", ylabel = \"T₂\", label = \"\")\nscatter!(1,1, markercolor = :green, markerstrokecolor = :green, marker = 2, label = \"converged\")\nscatter!(1,1, markercolor = :blue, markerstrokecolor = :blue, marker = 2, label = \"converged but hit bounds\")\nscatter!(1,1, markercolor = :red, markerstrokecolor = :red, marker = 2, label = \"not converged\")\n\nplt_2 = plot(plt_21, plt_22, layout = (1,2))\nplot(plt_1, plt_2, layout = grid(2,1, heights = [0.5, 0.5]), size=(700, 600))","category":"page"},{"location":"index.html#Optimal-control-problem","page":"Introduction","title":"Optimal control problem","text":"","category":"section"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"We consider the following optimal control problem ","category":"page"},{"location":"index.html","page":"Introduction","title":"Introduction","text":" left beginarrayll\n displaystyle min_xu int_t_0^t_f x(t) mathrm dt 1em\n textscdot x(t) = u(t) tin t_0 t_fmathrmae 05em\n phantommathrmsc u(t) in -11 tin t_0 t_f 05em\n phantommathrmsc x(t_0) = x_0 quad x(t_f) = x_f\n endarray right","category":"page"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"with x_0, t_0, x_f and t_f fixed. This problem is simple, and can be analytically solve without the use of numerical method. However, the goal is to solve this problem by indirect shooting. ","category":"page"},{"location":"index.html#Dependencies","page":"Introduction","title":"Dependencies","text":"","category":"section"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"All the numerical simulations to generate this documentation are performed with the following packages.","category":"page"},{"location":"index.html","page":"Introduction","title":"Introduction","text":"using Pkg\nPkg.status()","category":"page"},{"location":"2D-example.html#Indirect-method","page":"Indirect shooting","title":"Indirect method","text":"","category":"section"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"We introduce the pseudo-Hamiltonian ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":" h(xpp^0u) = p^0 x + p u","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"For the sake of simplicity, we consider in this notebook only the normal case, and we fix p^0 = -1. According to the Pontryagin maximum principle, the maximizing control is given by u(xp) to mathrmsign(p). This function is non-differentiable, and may lead to numerical issues. ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"using OptimalControl # hide\nusing Plots # hide\nusing ForwardDiff # hide\nusing DifferentialEquations # hide\nusing MINPACK # hide\n\nt0 = 0 # initial time\nx0 = 0 # initial state\ntf = 5 # final time\nxf = 0 # final state\n\n@def ocp begin # problem definition\n\n t ∈ [ t0, tf ], time\n x ∈ R, state\n u ∈ R, control\n\n x(t0) == x0\n x(tf) == xf\n\n ẋ(t) == u(t) \n\n ∫( x(t) ) → min\n\nend\nnothing # hide","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"using MINPACK\nfunction fsolve(f, j, x; kwargs...)\n try\n MINPACK.fsolve(f, j, x; kwargs...)\n catch e\n println(\"Erreur using MINPACK\")\n println(e)\n println(\"hybrj not supported. Replaced by hybrd even if it is not visible on the doc.\")\n MINPACK.fsolve(f, x; kwargs...)\n end\nend\nfunction fsolve(f, x; kwargs...)\n MINPACK.fsolve(f, x; kwargs...)\nend","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"Thanks to the control-toolbox, the flow varphi of the (true) Hamiltonian","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":" H(xp) = h(xp-1 u(xp)) = p^0 x + lvert p rvert ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"is given by the function textttFlow. The shooting function S colon mathbbR to mathbbR is defined by","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":" S(p_0) = pi big( varphi(t_0 x_0 p_0 t_f) big) - x_f","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"where pi (xp) = x is the classical x-space projection.","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"ϕ = Flow(ocp, (x,p) -> sign(p)) # flow with maximizing control \nπ((x,p)) = x; # projection on state space\n\nS(p0) = π( ϕ(t0, x0, p0, tf) ) - xf; # shooting function\nnle = p0 -> [S(p0[1])] # intermediate function\n\n# Plot\nplot(range(-7, 2, 500), S, xlim = [-7, 2])\nplot!([-7,2], [0,0], color = :black)\nplot!(xlabel = \"p0\", ylabel = \"S\", legend=false)","category":"page"},{"location":"2D-example.html#Finite-difference-method","page":"Indirect shooting","title":"Finite difference method","text":"","category":"section"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"The main goal now is to find the zero of S. To this purpose, we use the numerical solver texttthybrd1 given in the package textttMINPACKjl. If we don't provide the Jacobian J_S of S to the solver, the finite difference method is used to approximate it. ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"ξ = [-1.0] # initial guess\nS!(s, ξ) = (s[:] .= S(ξ[1]); nothing) # intermediate function\np0_sol = fsolve(S!, ξ, show_trace = true) # solve\nprintln(p0_sol)","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"sol = ϕ((t0, tf), x0, p0_sol.x) # get the optimal trajectory\nplot(sol) # plot","category":"page"},{"location":"2D-example.html#Automatic-differentiation-(wrong-way)","page":"Indirect shooting","title":"Automatic differentiation (wrong way)","text":"","category":"section"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"Now, we want to provide J_S to the solver, thanks to the textttForwardDiffjl package. This Jacobian is computed with the variational equation, and leads to a false result in our case.","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"Details. ( To be removed )","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"Denoting z = (xp), we have ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":" varphi(t_0 z_0 t_f) = z_0 + int_t_0^t_f vec Hbig(varphi(t_0 z_0 t)big) mathrm dt ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"If we assume that z_0 to varphi(t_0 z_0 t_f) is differentiable, we have ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":" fracpartial varphipartial z_0(t_0 z_0 t_f)cdot delta z_0 = delta z_0 + int_t_0^t_f vec Hbig(varphi(t_0 z_0 t)big)cdot left( fracpartial varphipartial z_0(t_0 z_0 t) cdot delta z_0 right) mathrm dt ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"and so, z_0 to fracpartial varphipartial z_0(t_0 z_0 t_f)cdot delta z_0 is solution of the variational equations","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":" fracpartial delta zpartial t(t) = vec Hbig(varphi(t_0 z_0 t_f)big) cdot delta z(t) qquad delta z(t_0) = delta z_0","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"In the studied optimal control problem, we have ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":" vec H(xp) = (mathrmsign(p) -1) ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"and so, we have vec H(z) = 0_2 almost everywhere, which implies","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":" fracpartial varphipartial z_0(t_0 z_0 t_f) cdot delta z_0 = mathrmexpbig((t_f-t_0) 0_2 big)cdot delta z_0 = delta z_0","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"The Jacobian of the shooting function is then given by ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":" S(p_0) = pi left( fracpartial varphipartial p_0(t_0 x_0 p_0 t_f) right) = pi left( fracpartial varphipartial z_0(t_0 z_0 t_f) cdot (01) right) = pi(01) = 0 ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"
\n
\n \n Details.\n
\n
\n

Denoting z=(x,p)z = (x,p)z=(x,p), we have

\n

φ(t0,z0,tf)=z0+∫t0tfH⃗(φ(t0,z0,t)) dt. \\varphi(t_0, z_0, t_f) = z_0 + \\int_{t_0}^{t_f} \\vec H\\big(\\varphi(t_0, z_0, t)\\big) \\,\\mathrm dt. φ(t0​,z0​,tf​)=z0​+∫t0​tf​​H(φ(t0​,z0​,t))dt.

\n

If we assume that z0→φ(t0,z0,tf)z_0 \\to \\varphi(t_0, z_0, t_f)z0​→φ(t0​,z0​,tf​) is differentiable, we have

\n

∂φ∂z0(t0,z0,tf)⋅δz0=δz0+∫t0tfH⃗′(φ(t0,z0,t))⋅(∂φ∂z0(t0,z0,t)⋅δz0) dt, \\frac{\\partial \\varphi}{\\partial z_0}(t_0, z_0, t_f)\\cdot \\delta z_0 = \\delta z_0 + \\int_{t_0}^{t_f} \\vec H'\\big(\\varphi(t_0, z_0, t)\\big)\\cdot \\left( \\frac{\\partial \\varphi}{\\partial z_0}(t_0, z_0, t) \\cdot \\delta z_0 \\right) \\,\\mathrm dt, ∂z0​∂φ​(t0​,z0​,tf​)⋅δz0​=δz0​+∫t0​tf​​H′(φ(t0​,z0​,t))⋅(∂z0​∂φ​(t0​,z0​,t)⋅δz0​)dt,

\n

and so, z0→∂φ∂z0(t0,z0,tf)⋅δz0z_0 \\to \\frac{\\partial \\varphi}{\\partial z_0}(t_0, z_0, t_f)\\cdot \\delta z_0z0​→∂z0​∂φ​(t0​,z0​,tf​)⋅δz0​ is solution of the variational equations

\n

∂δz∂t(t)=H⃗′(φ(t0,z0,tf))⋅δz(t),δz(t0)=δz0. \\frac{\\partial \\delta z}{\\partial t}(t) = \\vec H'\\big(\\varphi(t_0, z_0, t_f)\\big) \\cdot \\delta z(t), \\qquad \\delta z(t_0) = \\delta z_0.∂t∂δz​(t)=H′(φ(t0​,z0​,tf​))⋅δz(t),δz(t0​)=δz0​.

\n

In the studied optimal control problem, we have

\n

H⃗(x,p)=(sign(p),−1) \\vec H(x,p) = (\\mathrm{sign}(p), -1) H(x,p)=(sign(p),−1)

\n

and so, we have H⃗′(z)=02\\vec H'(z) = 0_2H′(z)=02​ almost everywhere, which implies

\n

∂φ∂z0(t0,z0,tf)⋅δz0=exp((tf−t0)02)⋅δz0=δz0. \\frac{\\partial \\varphi}{\\partial z_0}(t_0, z_0, t_f) \\cdot \\delta z_0 = \\mathrm{exp}\\big((t_f-t_0) 0_2 \\big)\\cdot \\delta z_0 = \\delta z_0.∂z0​∂φ​(t0​,z0​,tf​)⋅δz0​=exp((tf​−t0​)02​)⋅δz0​=δz0​.

\n

The Jacobian of the shooting function is then given by

\n

S′(p0)=π(∂φ∂p0(t0,x0,p0,tf))=π(∂φ∂z0(t0,z0,tf)⋅(0,1))=π(0,1)=0. S'(p_0) = \\pi \\left( \\frac{\\partial \\varphi}{\\partial p_0}(t_0, x_0, p_0, t_f) \\right) = \\pi \\left( \\frac{\\partial \\varphi}{\\partial z_0}(t_0, z_0, t_f) \\cdot (0,1) \\right) = \\pi(0,1) = 0. S′(p0​)=π(∂p0​∂φ​(t0​,x0​,p0​,tf​))=π(∂z0​∂φ​(t0​,z0​,tf​)⋅(0,1))=π(0,1)=0.

\n
\n
","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"ξ = [-1.0] # initial guess\nJS(ξ) = ForwardDiff.jacobian(p -> [S(p[1])], ξ) # compute jacobian by forward differentiation\nprintln(\"ξ = \", ξ[1])\nprintln(\"JS(ξ) : \", JS(ξ)[1])","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"However, the solver texttthybrd1 uses rank 1 approximation to actualize the Jacobian insted of compute it at each iteration, which imply that it still converges to the solution even if the given Jacobian is completely false.","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"JS!(js, ξ) = (js[:] .= JS(ξ); nothing) # intermediate function\np0_sol = fsolve(S!, JS!, ξ, show_trace = true) # solve\nprintln(p0_sol)","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"sol = ϕ((t0, tf), x0, p0_sol.x) # get the optimal trajectory\nplt = plot(sol) # plot","category":"page"},{"location":"2D-example.html#Automatic-differentiation-(good-way)","page":"Indirect shooting","title":"Automatic differentiation (good way)","text":"","category":"section"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"The goal is to provide the true Jacobian of S by using the textttForwardDiff package, and so we need to indicate to the solver that the dynamic of the system change when p = 0.","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"To understang why we need to give this information to the solver, see the following details. ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"Details. (To be removed)","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"The problem is that the Hamiltonian H is not differentiable everywhere due to the maximizing control. This control is bang-bang (u = 1 and u = -1). ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"Let now construct the two smooth Hamiltonians associated to these two controls ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":" H^+(xp) = h(xp-11) = -x + p qquad textand qquad H^-(xp) = h(xp-1-1) = -x - p ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"Their associated vector fields are given by ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":" vec H^+(xp) = (11) qquad textand qquad vec H^-(xp) = (-1 1) ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"and their associated flow correspond to ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":" varphi^+(t_0 z_0 t_f) = z_0 + left( beginarrayc 1 1 endarray right) (t_f -t_0)\n qquad textand qquad \n varphi^-(t_0 z_0 t_f) = z_0 + left( beginarrayc -1 phantom- 1 endarray right) (t_f -t_0)","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"If we assume that the optimal structure of the problem is negative then positive bangs, then the associated flow is defined by ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":" varphi(t_0 z_0 t_f) = varphi^+ big( t_1(z_0) varphi^-big(t_0 z_0 t_1(z_0)big) t_f big)","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"with the following condition ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":" pi_p big( varphi^-(t_0 z_0 t_1(z_0)) big) = 0","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"where pi_p(xp) = p is the classical p-space projection. By devlopping this last condition, an explicit form of the function t_1(cdot) is given by ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":" t_1(x_0 p_0) = t_0 - p_0","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"Finally, we have ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":" beginalign*\n fracpartial varphipartial z_0 \n = fracpartial varphi^+partial t_0 fracpartial t_1partial z_0 + fracvarphi^+partial z_0 left( fracpartial varphi^-partial z_0 + fracpartial varphi^-partial t_f fracpartial t_1partial z_0 right) \n = left( beginarrayc -1 -1 endarray right) left( beginarraycc0 -1 endarray right) + left( beginarraycc 1 0 0 1 endarray right) left\n left( beginarraycc 1 0 0 1 endarray right) + left( beginarrayc -1 phantom - 1 endarray right) left( beginarraycc0 -1 endarray right)\n right \n = left( beginarraycc 0 1 0 1 endarray right) + left( beginarraycc 1 0 0 1 endarray right) + left( beginarraycc 0 phantom -1 0 -1 endarray right) \n = left( beginarraycc 1 2 0 1 endarray right)\n endalign*","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"and so, we have that ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":" S(p_0) = pi left( fracpartial varphipartial p_0(t_0 x_0 p_0 t_f) right) = pi left( fracpartial varphipartial z_0(t_0 z_0 t_f) cdot (01) right) = pi(21) = 2","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"
\n
\n \n Details.\n
\n
\n

The problem is that the Hamiltonian HHH is not differentiable everywhere due to the maximizing control. This control is bang-bang (u=1u = 1u=1 and u=−1u = -1u=−1).

\n

Let now construct the two smooth Hamiltonians associated to these two controls

\n

H+(x,p)=h(x,p,−1,1)=−x+pandH−(x,p)=h(x,p,−1,−1)=−x−p. H^+(x,p) = h(x,p,-1,1) = -x + p \\qquad \\text{and} \\qquad H^-(x,p) = h(x,p,-1,-1) = -x - p. H+(x,p)=h(x,p,−1,1)=−x+pandH−(x,p)=h(x,p,−1,−1)=−x−p.

\n

Their associated vector fields are given by

\n

H⃗+(x,p)=(1,1)andH⃗−(x,p)=(−1,1), \\vec H^+(x,p) = (1,1) \\qquad \\text{and} \\qquad \\vec H^-(x,p) = (-1, 1), H+(x,p)=(1,1)andH−(x,p)=(−1,1),

\n

and their associated flow correspond to

\n

φ+(t0,z0,tf)=z0+(11)(tf−t0)andφ−(t0,z0,tf)=z0+(−1−1)(tf−t0). \\varphi^+(t_0, z_0, t_f) = z_0 + \\left( \\begin{array}{c} 1 \\\\ 1 \\end{array} \\right) (t_f -t_0)\n \\qquad \\text{and} \\qquad \n \\varphi^-(t_0, z_0, t_f) = z_0 + \\left( \\begin{array}{c} -1 \\\\ \\phantom{-} 1 \\end{array} \\right) (t_f -t_0).φ+(t0​,z0​,tf​)=z0​+(11​)(tf​−t0​)andφ−(t0​,z0​,tf​)=z0​+(−1−1​)(tf​−t0​).

\n

If we assume that the optimal structure of the problem is negative then positive bangs, then the associated flow is defined by

\n

φ(t0,z0,tf)=φ+(t1(z0),φ−(t0,z0,t1(z0)),tf), \\varphi(t_0, z_0, t_f) = \\varphi^+ \\big( t_1(z_0), \\varphi^-\\big(t_0, z_0, t_1(z_0)\\big), t_f \\big),φ(t0​,z0​,tf​)=φ+(t1​(z0​),φ−(t0​,z0​,t1​(z0​)),tf​),

\n

with the following condition

\n

πp(φ−(t0,z0,t1(z0)))=0, \\pi_p \\big( \\varphi^-(t_0, z_0, t_1(z_0)) \\big) = 0,πp​(φ−(t0​,z0​,t1​(z0​)))=0,

\n

where πp(x,p)=p\\pi_p(x,p) = pπp​(x,p)=p is the classical ppp-space projection. By devlopping this last condition, an explicit form of the function t1(⋅)t_1(\\cdot)t1​(⋅) is given by

\n

t1(x0,p0)=t0−p0. t_1(x_0, p_0) = t_0 - p_0.t1​(x0​,p0​)=t0​−p0​.

\n

Finally, we have

\n

∂φ∂z0=∂φ+∂t0∂t1∂z0+φ+∂z0(∂φ−∂z0+∂φ−∂tf∂t1∂z0)=(−1−1)(0−1)+(1001)[(1001)+(−1−1)(0−1)]=(0101)+(1001)+(0−10−1)=(1201) \\begin{align*}\n \\frac{\\partial \\varphi}{\\partial z_0} \n &= \\frac{\\partial \\varphi^+}{\\partial t_0} \\frac{\\partial t_1}{\\partial z_0} + \\frac{\\varphi^+}{\\partial z_0} \\left( \\frac{\\partial \\varphi^-}{\\partial z_0} + \\frac{\\partial \\varphi^-}{\\partial t_f} \\frac{\\partial t_1}{\\partial z_0} \\right) \\\\\n &= \\left( \\begin{array}{c} -1 \\\\ -1 \\end{array} \\right) \\left( \\begin{array}{cc}0 & -1 \\end{array} \\right) + \\left( \\begin{array}{cc} 1 & 0 \\\\ 0 & 1 \\end{array} \\right) \\left[\n \\left( \\begin{array}{cc} 1 & 0 \\\\ 0 & 1 \\end{array} \\right) + \\left( \\begin{array}{c} -1 \\\\ \\phantom - 1 \\end{array} \\right) \\left( \\begin{array}{cc}0 & -1 \\end{array} \\right)\n \\right] \\\\\n &= \\left( \\begin{array}{cc} 0 & 1 \\\\ 0 & 1 \\end{array} \\right) + \\left( \\begin{array}{cc} 1 & 0 \\\\ 0 & 1 \\end{array} \\right) + \\left( \\begin{array}{cc} 0 & \\phantom -1 \\\\ 0 & -1 \\end{array} \\right) \\\\ \n &= \\left( \\begin{array}{cc} 1 & 2 \\\\ 0 & 1 \\end{array} \\right)\n \\end{align*}∂z0​∂φ​​=∂t0​∂φ+​∂z0​∂t1​​+∂z0​φ+​(∂z0​∂φ−​+∂tf​∂φ−​∂z0​∂t1​​)=(−1−1​)(0​−1​)+(10​01​)[(10​01​)+(−1−1​)(0​−1​)]=(00​11​)+(10​01​)+(00​−1−1​)=(10​21​)​

\n

and so, we have that

\n

S′(p0)=π(∂φ∂p0(t0,x0,p0,tf))=π(∂φ∂z0(t0,z0,tf)⋅(0,1))=π(2,1)=2. S'(p_0) = \\pi \\left( \\frac{\\partial \\varphi}{\\partial p_0}(t_0, x_0, p_0, t_f) \\right) = \\pi \\left( \\frac{\\partial \\varphi}{\\partial z_0}(t_0, z_0, t_f) \\cdot (0,1) \\right) = \\pi(2,1) = 2.S′(p0​)=π(∂p0​∂φ​(t0​,x0​,p0​,tf​))=π(∂z0​∂φ​(t0​,z0​,tf​)⋅(0,1))=π(2,1)=2.

\n
\n
","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"To provide this change of dynamic to the solver, we need to use a callback during the integration that will execute the function textttaffect when textttcondition(xp) = 0.","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"For us, the condition is given by (xp) to p. For the textttaffect function, we use a global parameter alpha. This parameter will be set to pm 1 at the beginning of the integration and it sign will change with the textttaffect function. ","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"Thanks to the textttcontrol-toolbox package, the created callback can be easily pass to the integrator throught the textttFlow function.","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"global α # parameter: ̇p(t) = α with α = ±1\n\nfunction condition(z, t, integrator) # event when condition(x,p) == 0\n x,p = z\n return p\nend\n\nfunction affect!(integrator) # action when condition == 0 \n global α = -α\n nothing\nend\n\ncb = ContinuousCallback(condition, affect!) # callback \n\nφ_ = Flow(ocp, (x,p) -> α, callback = cb) # intermediate flow\n\nfunction φ(t0, x0, p0, tf; kwargs...) # flow\n global α = sign(p0)\n return φ_(t0, x0, p0, tf; kwargs...)\nend\n\nfunction φ((t0, tf), x0, p0; kwargs...) # flow for plot\n global α = sign(p0)\n return φ_((t0, tf), x0, p0; kwargs...)\nend\n\nShoot(p0) = π( φ(t0, x0, p0, tf) ) - xf # shooting function\n\nξ = [-1.0] # initial guess\nJShoot(ξ) = ForwardDiff.jacobian(p -> [Shoot(p[1])], ξ) # compute jacobian by forward differentiation\nprintln(\"ξ = \", ξ[1])\nprintln(\"JS(ξ) : \", JShoot(ξ)[1])\nShoot!(shoot, ξ) = (shoot[:] .= Shoot(ξ[1]); nothing) # intermediate function\nJShoot!(jshoot, ξ) = (jshoot[:] .= JShoot(ξ); nothing) # intermediate function\n\np0_sol = fsolve(Shoot!, JShoot!, ξ, show_trace = true) # solve\nprintln(p0_sol)","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"# get optimal trajectory\nsol_ = φ((t0, tf), x0, p0_sol.x[1], saveat=range(t0, tf, 500)) \n\n# plot\nsol = OptimalControl.OptimalControlSolution(sol_)\nt = sol.times\nx = sol.state\np = sol.costate\nu = sign ∘ p\n\nplt_x = plot(t, x, label = \"x\")\nplt_p = plot(t, p, label = \"p\")\nplt_u = plot(t, u, label = \"u\")\n\nplt_xp = plot(plt_x, plt_p, layout=(1, 2))\nplot(plt_xp, plt_u, layout = (2, 1))","category":"page"}] }