From 7802a4c5208ad0bccf871461325fe389c37c0528 Mon Sep 17 00:00:00 2001 From: "Documenter.jl" Date: Tue, 13 Aug 2024 07:21:14 +0000 Subject: [PATCH] build based on 1f3afa5 --- dev/.documenter-siteinfo.json | 2 +- ...e-8c496d8e.svg => 2D-example-1d185769.svg} | 192 +++++++++--------- ...e-7b1318a7.svg => 2D-example-367854fd.svg} | 192 +++++++++--------- ...e-126fd4a0.svg => 2D-example-66b689d0.svg} | 192 +++++++++--------- ...e-54bb81a9.svg => 2D-example-d58d3110.svg} | 58 +++--- dev/2D-example.html | 101 +++++++-- dev/index.html | 2 +- dev/search_index.js | 2 +- 8 files changed, 408 insertions(+), 333 deletions(-) rename dev/{2D-example-8c496d8e.svg => 2D-example-1d185769.svg} (80%) rename dev/{2D-example-7b1318a7.svg => 2D-example-367854fd.svg} (80%) rename dev/{2D-example-126fd4a0.svg => 2D-example-66b689d0.svg} (86%) rename dev/{2D-example-54bb81a9.svg => 2D-example-d58d3110.svg} (87%) diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json index 6732371..1575c64 100644 --- a/dev/.documenter-siteinfo.json +++ b/dev/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.10.4","generation_timestamp":"2024-08-12T16:51:38","documenter_version":"1.5.0"}} \ No newline at end of file +{"documenter":{"julia_version":"1.10.4","generation_timestamp":"2024-08-13T07:21:11","documenter_version":"1.5.0"}} \ No newline at end of file diff --git a/dev/2D-example-8c496d8e.svg b/dev/2D-example-1d185769.svg similarity index 80% rename from dev/2D-example-8c496d8e.svg rename to dev/2D-example-1d185769.svg index efda61a..5945eb2 100644 --- a/dev/2D-example-8c496d8e.svg +++ b/dev/2D-example-1d185769.svg @@ -1,118 +1,118 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/2D-example-7b1318a7.svg b/dev/2D-example-367854fd.svg similarity index 80% rename from dev/2D-example-7b1318a7.svg rename to dev/2D-example-367854fd.svg index cd4051d..a2d0cab 100644 --- a/dev/2D-example-7b1318a7.svg +++ b/dev/2D-example-367854fd.svg @@ -1,118 +1,118 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/2D-example-126fd4a0.svg b/dev/2D-example-66b689d0.svg similarity index 86% rename from dev/2D-example-126fd4a0.svg rename to dev/2D-example-66b689d0.svg index b65c004..06c875f 100644 --- a/dev/2D-example-126fd4a0.svg +++ b/dev/2D-example-66b689d0.svg @@ -1,118 +1,118 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/2D-example-54bb81a9.svg b/dev/2D-example-d58d3110.svg similarity index 87% rename from dev/2D-example-54bb81a9.svg rename to dev/2D-example-d58d3110.svg index 52cec58..98e6f33 100644 --- a/dev/2D-example-54bb81a9.svg +++ b/dev/2D-example-d58d3110.svg @@ -1,44 +1,44 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/2D-example.html b/dev/2D-example.html index 644eccd..4fba86f 100644 --- a/dev/2D-example.html +++ b/dev/2D-example.html @@ -26,19 +26,56 @@ # 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
Results of Nonlinear Solver Algorithm
+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.165962
+     2     9.318967e-01     1.069263e+00         0.240037
+     3     3.908179e-08     2.171079e-01         0.001056
+     4     3.126006e-08     3.818465e-16         0.000929
+     5     2.602624e-08     7.541206e-17         0.000943
+     6     1.122729e-07     1.864781e-15         0.001027
+     7     1.251064e-08     1.228963e-15         0.000997
+     8     1.071865e-08     6.960129e-18         0.000961
+     9     3.372047e-08     2.490162e-16         0.000935
+    10     6.149611e-09     1.433785e-16         0.000994
+    11     2.836331e-08     2.624355e-17         0.000956
+    12     2.550366e-08     1.772446e-17         0.000943
+    13     2.794534e-08     1.292187e-17         0.001848
+    14     2.543587e-08     1.364932e-17         0.000925
+    15     6.042093e-09     4.286522e-19         0.000988
+    16     3.036848e-08     7.912620e-17         0.000900
+    17     2.599380e-08     5.504418e-17         0.000942
+    18     2.518026e-08     1.434518e-18         0.000932
+    19     6.005500e-09     5.041160e-20         0.000986
+    20     3.036868e-08     7.817582e-17         0.000911
+    21     2.601931e-08     5.449264e-17         0.000938
+    22     2.521370e-08     1.406704e-18         0.000917
+    23     5.969729e-09     4.888090e-20         0.000977
+    24     3.036827e-08     7.723669e-17         0.000916
+    25     2.604422e-08     5.394382e-17         0.000920
+    26     2.524639e-08     1.379606e-18         0.000928
+    27     5.934758e-09     4.741150e-20         0.000997
+    28     3.036860e-08     7.634290e-17         0.000907
+    29     2.606871e-08     5.342255e-17         0.000920
+    30     2.527841e-08     1.353711e-18         0.000917
+    31     5.900549e-09     4.601691e-20         0.000984
+    32     3.036849e-08     7.546245e-17         0.000914
+    33     2.609265e-08     5.290603e-17         0.000929
+    34     2.530975e-08     1.328510e-18         0.000915
+Results of Nonlinear Solver Algorithm
  * Algorithm: Modified Powell
  * Starting Point: [-1.0]
  * Zero: [-2.500000007612481]
  * 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.450148 seconds
+ * Total time: 0.437277 seconds
  * Function Calls: 34
- * Jacobian Calls (df/dx): 2

!!! Don't show the trace !!! 

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. \]

+ * 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. \]

Details. @@ -92,17 +129,31 @@ println("ξ = ", ξ[1]) println("JS(ξ) : ", JS(ξ)[1])
ξ = -1.0
 JS(ξ) : 0.0

However, the solver $\texttt{hybrd1}$ 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.

JS!(js, ξ) = (js[:] .= JS(ξ); nothing)                      # intermediate function
-p0_sol = fsolve(S!, JS!, ξ, show_trace = true)              # solve
Results of Nonlinear Solver Algorithm
+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.006238
+     2     5.000000e+00     1.000000e+04         0.060365
+     3     5.000000e+00     3.906250e+03         0.000146
+     4     5.000000e+00     3.906250e+03         0.001357
+     5     5.000000e+00     3.906250e+03         0.000099
+     6     5.000000e+00     5.493164e+02         0.000080
+     7     5.000000e+00     7.724762e+01         0.000100
+     8     9.550781e-01     1.086295e+01         0.001118
+     9     2.633355e-08     2.280435e-01         0.000954
+    10     4.034347e-08     1.733639e-16         0.000994
+    11     7.549755e-09     6.346772e-17         0.000989
+Results of Nonlinear Solver Algorithm
  * Algorithm: Modified Powell (User Jac, Expert)
  * Starting Point: [-1.0]
  * Zero: [-2.500000011044682]
  * 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.106551 seconds
+ * Total time: 0.072453 seconds
  * Function Calls: 11
- * Jacobian Calls (df/dx): 2

!!! Don't show the trace !!!

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) + * 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) \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} @@ -193,16 +244,40 @@ Shoot!(shoot, ξ) = (shoot[:] .= Shoot(ξ[1]); nothing) # intermediate function JShoot!(jshoot, ξ) = (jshoot[:] .= JShoot(ξ); nothing) # intermediate function -p0_sol = fsolve(Shoot!, JShoot!, ξ, show_trace = true) # solve

Results of Nonlinear Solver Algorithm
+p0_sol = fsolve(Shoot!, JShoot!, ξ, show_trace = true)      # solve
+println(p0_sol)
ξ = -1.0
+JS(ξ) : 2.0000000000000018
+Iter     f(x) inf-norm    Step 2-norm      Step time
+------   --------------   --------------   --------------
+     1     3.000000e+00     0.000000e+00         1.298577
+     2     6.994405e-15     2.250000e+00         0.007679
+     3     4.551914e-15     1.262177e-29         0.000190
+     4     3.441691e-15     1.774937e-30         0.000148
+     5     2.775558e-15     1.972152e-31         0.000149
+     6     3.441691e-15     1.972152e-31         0.000134
+     7     2.775558e-15     1.972152e-31         0.000142
+     8     2.997602e-15     1.774937e-30         0.000297
+     9     1.443290e-15     1.972152e-31         0.000147
+    10     2.997602e-15     7.888609e-31         0.000139
+    11     2.997602e-15     1.972152e-31         0.000145
+    12     2.997602e-15     1.972152e-31         0.000275
+    13     2.997602e-15     1.972152e-31         0.000147
+    14     1.443290e-15     1.972152e-31         0.000130
+    15     2.343750e-02     1.373291e-04         0.000140
+    16     2.997602e-15     1.373291e-04         0.000131
+    17     2.997602e-15     1.972152e-31         0.000140
+    18     1.443290e-15     1.972152e-31         0.000132
+    19     1.464844e-03     5.364418e-07         0.000140
+Results of Nonlinear Solver Algorithm
  * Algorithm: Modified Powell (User Jac, Expert)
  * Starting Point: [-1.0]
  * Zero: [-2.5]
  * 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.347241 seconds
+ * Total time: 1.308996 seconds
  * Function Calls: 19
- * Jacobian Calls (df/dx): 3

!!! Don't show the trace !!!

# get optimal trajectory
+ * Jacobian Calls (df/dx): 3
# get optimal trajectory
 sol_ = φ((t0, tf), x0, p0_sol.x[1], saveat=range(t0, tf, 500))
 
 # plot
@@ -217,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/index.html b/dev/index.html index fa1b662..fa6dde3 100644 --- a/dev/index.html +++ b/dev/index.html @@ -4,4 +4,4 @@ \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.

+ \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.

diff --git a/dev/search_index.js b/dev/search_index.js index c6410ac..0e9151b 100644 --- a/dev/search_index.js +++ b/dev/search_index.js @@ -1,3 +1,3 @@ var documenterSearchIndex = {"docs": -[{"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":"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":"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","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"!!! Don't show the trace !!! ","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","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"!!! Don't show the trace !!!","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","category":"page"},{"location":"2D-example.html","page":"Indirect shooting","title":"Indirect shooting","text":"!!! Don't show the trace !!!","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":"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":"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"}] }