From 75f14136a6716ce93d441d9170e0f39fac3ddd7c Mon Sep 17 00:00:00 2001 From: "Documenter.jl" Date: Mon, 12 Aug 2024 16:51:42 +0000 Subject: [PATCH] build based on 98a9563 --- dev/.documenter-siteinfo.json | 2 +- dev/2D-example-126fd4a0.svg | 118 ++++++++++++++++++ dev/2D-example-54bb81a9.svg | 44 +++++++ dev/2D-example-7b1318a7.svg | 118 ++++++++++++++++++ dev/2D-example-8c496d8e.svg | 118 ++++++++++++++++++ dev/2D-example.html | 225 +++++++++++++++++++++++++++++++--- dev/index.html | 7 +- dev/objects.inv | Bin 308 -> 353 bytes dev/search_index.js | 2 +- 9 files changed, 613 insertions(+), 21 deletions(-) create mode 100644 dev/2D-example-126fd4a0.svg create mode 100644 dev/2D-example-54bb81a9.svg create mode 100644 dev/2D-example-7b1318a7.svg create mode 100644 dev/2D-example-8c496d8e.svg diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json index 00c8b81..6732371 100644 --- a/dev/.documenter-siteinfo.json +++ b/dev/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.10.4","generation_timestamp":"2024-08-12T14:03:47","documenter_version":"1.5.0"}} \ No newline at end of file +{"documenter":{"julia_version":"1.10.4","generation_timestamp":"2024-08-12T16:51:38","documenter_version":"1.5.0"}} \ No newline at end of file diff --git a/dev/2D-example-126fd4a0.svg b/dev/2D-example-126fd4a0.svg new file mode 100644 index 0000000..b65c004 --- /dev/null +++ b/dev/2D-example-126fd4a0.svg @@ -0,0 +1,118 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/2D-example-54bb81a9.svg b/dev/2D-example-54bb81a9.svg new file mode 100644 index 0000000..52cec58 --- /dev/null +++ b/dev/2D-example-54bb81a9.svg @@ -0,0 +1,44 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/2D-example-7b1318a7.svg b/dev/2D-example-7b1318a7.svg new file mode 100644 index 0000000..cd4051d --- /dev/null +++ b/dev/2D-example-7b1318a7.svg @@ -0,0 +1,118 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/2D-example-8c496d8e.svg b/dev/2D-example-8c496d8e.svg new file mode 100644 index 0000000..efda61a --- /dev/null +++ b/dev/2D-example-8c496d8e.svg @@ -0,0 +1,118 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/2D-example.html b/dev/2D-example.html index c0f152f..644eccd 100644 --- a/dev/2D-example.html +++ b/dev/2D-example.html @@ -1,21 +1,10 @@ -2D example · Geometric preconditioner
GitHub

Optimal control problem

We consider the following optimal control problem

\[ \left\{ \begin{array}{ll} - \displaystyle \min_{x,u} \int_{t_0}^{t_f} x(t) ~\mathrm dt \\[1em] - \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.

Indirect method

We introduce the pseudo-Hamiltonian

\[ h(x,p,p^0,u) = p^0 x + p u.\]

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(x,p) \to \mathrm{sign}(p)$. This function is non-differentiable, and may lead to numerical issues.

Let us import the necessary package and define the studied optimal control problem with some fixed initial and final time and state values.

using OptimalControl
-using Plots
-using ForwardDiff
-using DifferentialEquations
-using MINPACK
-
-t0 = 0
-x0 = 0
-tf = 5
-xf = 0                      # initial and final time and state
-
-@def ocp begin              # problem definition
+Indirect shooting · Geometric preconditioner
GitHub
Indirect method
We introduce the pseudo-Hamiltonian 
\[    h(x,p,p^0,u) = p^0 x + p u.\]
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(x,p) \to \mathrm{sign}(p)$. This function is non-differentiable, and may lead to numerical issues.  
t0 = 0                                                      # initial time
+x0 = 0                                                      # initial state
+tf = 5                                                      # final time
+xf = 0                                                      # final state
+
+@def ocp begin                                              # problem definition
 
     t ∈ [ t0, tf ], time
     x ∈ R, state
@@ -28,4 +17,204 @@
 
     ∫( x(t) ) → min
 
-end

+end

Thanks to the control-toolbox, the flow $\varphi$ of the (true) Hamiltonian

\[ H(x,p) = h(x,p,-1, u(x,p)) = p^0 x + \lvert p \rvert \]

is given by the function $\texttt{Flow}$. The shooting function $S \colon \mathbb{R} \to \mathbb{R}$ is defined by

\[ S(p_0) = \pi \big( \varphi(t_0, x_0, p_0, t_f) \big) - x_f\]

where $\pi (x,p) = x$ is the classical $x$-space projection.

ϕ = Flow(ocp, (x,p) -> sign(p))                             # flow with maximizing control
+π((x,p)) = x;                                               # projection on state space
+
+S(p0) = π( ϕ(t0, x0, p0, tf) ) - xf;                        # shooting function
+nle = p0 -> [S(p0[1])]                                      # intermediate function
+
+# 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
+S!(s, ξ) = (s[:] .= S(ξ[1]); nothing)                       # intermediate function
+p0_sol = fsolve(S!, ξ, show_trace = true)                   # solve
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
+ * 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. \]

+
+ + Details. +
+
+

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

+

φ(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.

+

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

+

∂φ∂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,

+

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_0 is solution of the variational equations

+

∂δ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.

+

In the studied optimal control problem, we have

+

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

+

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

+

∂φ∂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.

+

The Jacobian of the shooting function is then given by

+

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.

+
+
ξ = [-1.0]                                                  # initial guess
+JS(ξ) = ForwardDiff.jacobian(p -> [S(p[1])], ξ)             # compute jacobian by forward differentiation
+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
+ * 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
+ * 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) + \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} + &= \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) \\ + &= \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[ + \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) + \right] \\ + &= \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) \\ + &= \left( \begin{array}{cc} 1 & 2 \\ 0 & 1 \end{array} \right) + \end{align*}\]

and so, we have that

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

+
+ + Details. +
+
+

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

+

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

+

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.

+

Their associated vector fields are given by

+

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),

+

and their associated flow correspond to

+

φ+(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) + \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

+

φ(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),

+

with the following condition

+

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

+

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

+

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

+

Finally, we have

+

∂φ∂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*} + \frac{\partial \varphi}{\partial z_0} + &= \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) \\ + &= \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[ + \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) + \right] \\ + &= \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) \\ + &= \left( \begin{array}{cc} 1 & 2 \\ 0 & 1 \end{array} \right) + \end{align*}

+

and so, we have that

+

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.

+
+

To provide this change of dynamic to the solver, we need to use a callback during the integration that will execute the function $\texttt{affect!}$ when $\texttt{condition(x,p)} = 0$.

For us, the condition is given by $(x,p) \to p$. For the $\texttt{affect!}$ 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 $\texttt{affect!}$ function.

Thanks to the $\texttt{control-toolbox}$ package, the created callback can be easily pass to the integrator throught the $\texttt{Flow}$ function.

global α                                                    # parameter: ̇p(t) = α with α = ±1
+
+function condition(z, t, integrator)                        # event when condition(x,p) == 0
+    x,p = z
+    return p
+end
+
+function affect!(integrator)                                # action when condition == 0
+    global α = -α
+    nothing
+end
+
+cb = ContinuousCallback(condition, affect!)                 # callback
+
+φ_ = Flow(ocp, (x,p) -> α, callback = cb)                   # intermediate flow
+
+function φ(t0, x0, p0, tf; kwargs...)                       # flow
+    global α = sign(p0)
+    return φ_(t0, x0, p0, tf; kwargs...)
+end
+
+function φ((t0, tf), x0, p0; kwargs...)                     # flow for plot
+    global α = sign(p0)
+    return φ_((t0, tf), x0, p0; kwargs...)
+end
+
+Shoot(p0) =  π( φ(t0, x0, p0, tf) ) - xf                    # shooting function
+
+ξ = [-1.0]                                                  # initial guess
+JShoot(ξ) = ForwardDiff.jacobian(p -> [Shoot(p[1])], ξ)     # compute jacobian by forward differentiation
+println("ξ = ", ξ[1])
+println("JS(ξ) : ", JShoot(ξ)[1])
+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
+ * 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
+ * Function Calls: 19
+ * Jacobian Calls (df/dx): 3

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

# get optimal trajectory
+sol_ = φ((t0, tf), x0, p0_sol.x[1], saveat=range(t0, tf, 500))
+
+# plot
+sol = OptimalControl.OptimalControlSolution(sol_)
+t = sol.times
+x = sol.state
+p = sol.costate
+u = sign ∘ p
+
+plt_x = plot(t, x, label = "x")
+plt_p = plot(t, p, label = "p")
+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
diff --git a/dev/index.html b/dev/index.html index 050c6b8..fa1b662 100644 --- a/dev/index.html +++ b/dev/index.html @@ -1,2 +1,7 @@ -Introduction · Geometric preconditioner
GitHub
+Introduction · Geometric preconditioner
GitHub

Optimal control problem

We consider the following optimal control problem

\[ \left\{ \begin{array}{ll} + \displaystyle \min_{x,u} \int_{t_0}^{t_f} x(t) ~\mathrm dt \\[1em] + \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.

diff --git a/dev/objects.inv b/dev/objects.inv index 626948889883203a42a6c31cd7a6e8828cccc1e9..9f69a162699b4c96815c766c736ed12982e82209 100644 GIT binary patch delta 222 zcmV<403rXh0^tIXi+`n%NeTin5JmTWiUJjM5y4)mJ3$0jo*^b>cR`YBVu^ZtjXgMk z1Gp{zd;cdCES_L5c5c~#ngFH4xD@YT9VD_n7$FW>zEpV9zH<1OzF%9!QbZ%LFPVix zk!>m=ubncSY)(~?w+{VH;`A2Xb7`*+4FL@S4Uw%Vkf5@(+gD-4U{HCKlI^qVWzp4F zb8zRN=CXF|kDl&r4Zb)Cdo5(TNkTgPl^Px+he93odGuNOUNOk0on(;Dpi9C_66!KY YWuGP$%=X_!t>ep}6sz}K04jjkc>iK%c>n+a delta 177 zcmV;i08anm0<;2 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"}] }