Notebook source code: examples/irrigation_hybrid/hybrid.ipynb
Run the notebook yourself on binder Binder badge

Singular arc and non-smooth dynamics in optimal control

An example

  • Author: Olivier Cots

  • Date: February 2022


Abstract.

We consider an optimal control problem from irrigation application which has a non-smooth dynamics. The problem may be written in Mayer form with a control-system affine with respect to the control and with points of non differentiability with respect to the state. The peculiar feature is that the optimal solution contains a singular arc which appears at a non differentiable point (called corner point) of the dynamics. We present a direct shooting formulation to solve this optimal control problem when the structure is known a priori and present a regularization process combined with a continuation technique to recover the fact that the singular arc is located at a corner point.

Main goal

Retrieve, by regularization and continuation, the fact the singular arc of the optimal solution is located at a corner point of the dynamics.

Problem formulation

We consider the following optimal control problem in Mayer form, from crop irrigation optimization application (see [Boumaza et al. 2020]):

\[\begin{split}\left\{ \begin{array}{l} \min J(x(\cdot), u(\cdot)) \overset{\mathrm{def}}{=} -B(T) \\[1.0em] \dot{S}(t) = -k_1\big(\varphi(t)\, K_S(S(t)) + (1-\varphi(t))\, K_R(S(t)) \big) + k_1k_2 u(t), \\[0.5em] \dot{B}(t) = \varphi(t)\, K_S(S(t))\, f(B(t)), \\[0.5em] \dot{V}(t) = u(t), \quad u(t) \in [0, 1], \quad t \in [0, T] \text{ a.e.}, \\[1.0em] S(0) = 1, \quad B(0)=B_0 >0, \quad V(0)=0, \quad V(T) = \bar{V}, \end{array} \right.\end{split}\]

where the final time \(T\) (the harvesting date) is fixed and where \(\bar{V} = F_\mathrm{max} / \bar{Q}\) is given. The states \(S(t)\), \(B(t)\) and \(V(t)\) stand respectively for the relative soil humidity in the root zone, the crop biomass at time \(t\) and the normalized total water delivered to the crop during the time interval \([0,t]\). We define the state vector \(x \overset{\mathrm{def}}{=} (S, B, V)\). The control variable \(u(t) = F(t) / F_\mathrm{max}\) is the ratio of the input water flow rate \(F(t)\) at time \(t\) over the maximal flow \(F_\mathrm{max}\) that the irrigation allows. See [Boumaza et al. 2020] for a more detailed description of the problem.

Remark. The final constaint \(V(T) = \bar{V}\) may be seen as a \(L^1\) constraint on the control.

Assumption 1. The functions \(K_S\) and \(K_R\) are piecewise linear non decreasing from \([0,1]\) to \([0,1]\) and are given by the following expressions

\[\begin{split}K_S(S) \overset{\mathrm{def}}{=} \left\{ \begin{array}{cl} 0 & \text{if } S \in [0, S_w] \\[0.3em] \displaystyle \frac{S-S_w}{S^\star-S_w} & \text{if } S \in [S_w, S^\star] \\[0.3em] 1 & \text{if } S \in [S^\star, 1] \end{array} \right. \quad \text{and} \quad K_R(S) \overset{\mathrm{def}}{=} \left\{ \begin{array}{cl} 0 & \text{if } S \in [0, S_h] \\[0.3em] \displaystyle \frac{S-S_h}{1-S_h} & \text{if } S \in [S_h, 1] \end{array} \right.\end{split}\]

where \(0 < S_h < S_w < S^\star < 1\). We shall say that \(S\) is a corner point of \(K_S\), resp. \(K_R\) when the function is non differentiable at \(S\).

3d504f5adaea4bbdac65809a81535c46

Figure: Graphs of the functions \(K_S\) and \(K_R\).

Assumption 2. \(\varphi\) is a \(C^1\) increasing function from \([0,T]\) to \([0,1]\) with \(\varphi(0)=0\) and \(\varphi(T)=1\). We consider

\[\varphi(t) = \left(\frac{t}{T}\right)^\alpha.\]

Assumption 3. \(k_1\), \(k_2\) are positive parameters with \(k_2 \ge 1\).

Assumption 4. The function \(f\) is the logistic law

\[f(B) = r B \left( 1 - \frac{B}{B_\mathrm{max}} \right),\]

with \(B_\mathrm{max} > B_0\).

The condition \(k_2\ge 1\) is a controllability assumption, in the sense that it allows the variable \(S\) to stay equal to \(1\) with the constant control \(u = 1/k_2\). This is desired since the state \(S\) in naturally subject to the (state) constraint \(S(t) \le 1\), \(t \in [0,T]\). We shall consider the set of admissible controls \(\mathcal{U}\) as measurable functions \(u(\cdot)\) taking values in \([0,1]\) such that the solution of the control-system satisfies the state constraint.

The optimal control problem we have defined has several difficulties to tackle:

  1. the control appears linearly in the (extended) control-system so the optimal solution may be composed of bang arcs (such that \(u(t) \in \{0, 1\}\)) and singular arcs (such that \(u(t) \in (0,1)\), see [Bonnard and Chyba 2003]) with possibly discontinuities at the junctions of the bang and singular arcs. Hence, the optimal control is possibly non-continuous (but still piecewise smooth). See [Bonnard et al. 2015], [Bonnard and Cots 2014], [Bonnard et al. 2020] for the application of indirect shooting and homotopy techniques for the computation of optimal control laws with bang and singular arcs.

  2. the problem has a state constraint (of order 1) \(S(t) \in [0,1]\) and so the optimal solution may also contain boundary arcs. See [Cots 2017], [Cots et al. 2018] for the application of direct methods, indirect shooting and homotopy techniques for the computation of optimal control laws with bang and boundary arcs.

  3. the problem dynamics is non-differentiable when \(S \in \{S_w, S^\star, S_h\}\) and so we need a non-smooth version of the Pontryagin Maximum Principle [Clarke 2013] or a Maximum Principle for hybrid optimal control problems (with state constraints if the optimal solution has boundary arcs) [Caines et al. 2006, Dmitruk and Kaganovich 2008], [Haberkorn and Trélat 2011], [Sussmann 1999], [Sussmann 1999 (GOC)].

Remark. Note that according to [Boumaza et al. 2020], the singular arcs only occur at corner points.

Main difficulty

Due to difficulties 1 and 2, the optimal control may have a complex structure. In this work, we assume we know the structure, that is the sequence of bang, singular and boundary arcs. Note that there are no boundary arcs in the optimal solution so we can remove this constraint. The main difficulty that we have to deal with comes from the non differentiability of the dynamics at the point \(S^\star\) since the optimal solution has only one singular arc, which is located at \(S=S^\star\).

The state domain is \(M \overset{\mathrm{def}}{=} [0,1] \times [0, B_\mathrm{max}] \times \mathrm{R}_+\). One can stratify \(M\) according to

\[M = M^- \cup M^0 \cup M^+,\]

where the strata are defined by

\[\begin{split}\begin{aligned} M^- & \overset{\mathrm{def}}{=} \{ x = (S,B,V) \in M ~|~ S < S^\star \}, \\ M^0 & \overset{\mathrm{def}}{=} \{ x = (S,B,V) \in M ~|~ S = S^\star \}, \\ M^+ & \overset{\mathrm{def}}{=} \{ x = (S,B,V) \in M ~|~ S > S^\star \}. \end{aligned}\end{split}\]

Note that \(M^0\) is of dimension 2 while \(M^-\) and \(M^+\) are of dimension 3. The optimal control problem may be seen has an hybrid optimal control problem, with different dynamics on each stratum. The main difficulty comes from the fact that the optimal solution has a singular arc located on the stratum \(M^0\) which is at the interface between \(M^-\) and \(M^+\), two strata of higher dimension.

Numerical preliminaries

For the numercial experiments, we consider the following normalized values for the parameters.

\(T\)

\(k_1\)

\(k_2\)

\(S^\star\)

\(S_w\)

\(S_h\)

\(F_\mathrm{max}\)

\(\bar{Q}\)

\(\alpha\)

\(r\)

\(B_\mathrm{max}\)

\(B_0\)

\(1\)

\(2.1\)

\(5\)

\(0.7\)

\(0.4\)

\(0.2\)

\(1.2\)

\(0.1\)

\(4\)

\(25\)

\(1\)

\(0.0005\)

[1]:
# Packages
#
import nutopy as nt # For indirect shooting and flows
import numpy as np  # For numerics
from scipy.optimize import minimize, NonlinearConstraint, LinearConstraint # For direct shooting
#
# Personal files
import plottings    # for plots of the solutions
import wrappers     # to get the Hamiltonians and the singular control
import bdd          # for saving data
import continuation # a simple discrete continuation method
[2]:
# Parameters
T     = 1.0
k1    = 2.1
k2    = 5.0
Ss    = 0.7 # Sstar
Sw    = 0.4
Sh    = 0.2
Fmax  = 1.2
Qbar  = 0.1
alpha = 4.0
r     = 25.0
Bmax  = 1
B0    = 0.0005

t0    = 0.0
x0_   = np.array([1.0, B0, 0.0])
Vbar  = Qbar/Fmax
pzero = -1.0

# The different values of eta for the regularization of KS
eta_init = 0.3
eta_span = np.arange(eta_init, 0.0, -0.01)
eta_span = np.append(eta_span, 0.001)

# Initialize data
data_file = 'data.json'
restart   = False # restart or not the computations
data      = bdd.Data({'eta_span': eta_span.tolist()}, data_file, restart)

Back to top

Optimal solution via direct shooting

If the optimal solution is such that \(S(t)\) is constant for times \(t \in [a, b]\) with \(a < b\) then on \((a,b)\), we have \(\dot{S}(t) = 0\). If we denote by \(\bar{S}\) the constant, then, the optimal control \(u\) must satisfy \(u(t) = \xi(t, \bar{S})\) on \((a,b)\), where

\[\xi(t, S) \overset{\mathrm{def}}{=} \frac{1}{k_2} \big( \varphi(t) K_S(S) + (1-\varphi(t)) K_R(S)\big).\]

Let us introduce the control law \(u^\star(t) \overset{\mathrm{def}}{=} \xi(t, S^\star)\). Then, since \(K_S(S^\star) = 1\), we have

\[u^\star(t) = \frac{1}{k_2} \big( \varphi(t) + (1-\varphi(t)) K_R(S^\star)\big).\]

Goal 1. We assume that the optimal solution is of the form

\[\begin{split}u(t) = \left\{ \begin{array}{lll} 0 & \text{if} & t \in [0, t_1], \\[0.5em] 1 & \text{if} & t \in (t_1, t_2], \\[0.5em] u^\star(t) & \text{if} & t \in (t_2, t_3], \\[0.5em] 0 & \text{if} & t \in (t_3, T] \end{array} \right.\end{split}\]

with \(S(t)=S^\star\) for \(t \in [t_2,t_3]\) and we aim to compute numerically, by direct shooting [Bock and Plitt 1984], the optimal switching times \(t_1\), \(t_2\) and \(t_3.\)

Exponential map. We introduce the notation

\[\dot{x}(t) = F_0(t, x(t)) + u(t)\, F_1, \quad x = (S, B, V),\]

\(F_1\) being a constant vector field. To solve the optimal control problem by direct shooting, we define the following three different vector fields:

\[F^{\min}(t, x) \overset{\mathrm{def}}{=} F_0(t, x), \quad F^{\max}(t, x) \overset{\mathrm{def}}{=} F_0(t, x) + F_1 \quad \text{and} \quad F^\star(t, x) \overset{\mathrm{def}}{=} F_0(t, x) + u^\star(t)\, F_1.\]

We define also the exponential map \(\exp_F(t_1, t_0, x_0)\) as the solution at time \(t_1\) of \(\dot x(t) = F(t, x(t))\), \(x(t_0) = x_0\).

[3]:
# KR and KS: see Assumption 1
def KS(S):
    y = 0.
    if(S>Sw and S<Ss):
        y = (S-Sw)/(Ss-Sw)
    elif(S>=Ss):
        y = 1.
    return y

def KR(S):
    y = 0.
    if(S>Sh):
        y = (S-Sh)/(1.-Sh)
    return y

# phi: see Assumption 2
def phi(t):
    return (t/T)**alpha

# f: see Assumption 4
def f(B):
    return r*B*(1.0-B/Bmax)

# F0 and F1: see the optimal control problem dynamics
def F0(t, x):
    S = x[0]
    B = x[1]
    return np.array([k1*(-phi(t)*KS(S)-(1.0-phi(t))*KR(S)), phi(t)*KS(S)*f(B), 0.0])

F1 = np.array([k1*k2, 0., 1.])

# umin = 0
@nt.tools.vectorize(vvars=(1,))
def umin(t):
    return 0.0

# Fmin = F0
def Fmin(t, x):
    return F0(t, x)

# umax
@nt.tools.vectorize(vvars=(1,))
def umax(t):
    return 1.0

# Fmax = F0 + umax F1 = F0 + F1
def Fmax(t, x):
    return F0(t, x) + F1

# ustar(t) = xi(t, Ss)
@nt.tools.vectorize(vvars=(1,))
def ustar(t):
    return (phi(t)+(1.0-phi(t))*KR(Ss))/k2

# Fstar = F0 + ustar F1
def Fstar(t, x):
    return F0(t, x) + ustar(t) * F1

# Exponential map
def exponential(F, t1, t0, x0):
    sol = nt.ivp.exp(F, t1, t0, x0)
    return sol.xf

Direct shooting and optimization problem. We introduce the following endpoint mapping which assigns to the switching times \((t_1, t_2, t_3)\) the final state \(x=(S,B,V)\) at time \(T\), departing from \(x_0\) at time \(t_0\):

\[\Theta(t_1,t_2,t_3) \overset{\mathrm{def}}{=} \exp_{F^{\min}}(T, t_3, \cdot) \circ \exp_{F^\star}(t_3, t_2, \cdot) \circ \exp_{F^{\max}}(t_2, t_1, \cdot) \circ \exp_{F^{\min}}(t_1, t_0, \cdot)(x_0),\]

\(t_0\), \(x_0\) and \(T\) being given. We say that the trajectory has a structure of the form \(\sigma^{\min}\sigma^{\max}\sigma^{\star}\sigma^{\min}\). We define the following map which gives an intermediate point:

\[\Psi(t_1,t_2) \overset{\mathrm{def}}{=} \exp_{F^{\max}}(t_2, t_1, \cdot) \circ \exp_{F^{\min}}(t_1, t_0, \cdot)(x_0),\]

such that \(\Theta(t_1,t_2,t_3) = \exp_{F^{\min}}(T, t_3, \cdot) \circ \exp_{F^\star}(t_3, t_2, \cdot)(\Psi(t_1,t_2))\). We define also the following canonical projections: \(\pi_S(x) = S\), \(\pi_B(x) = B\) and \(\pi_V(x) = V\).

Finally, we introduce the corresponding optimization problem where the decision variable \(y \overset{\mathrm{def}}{=} (t_1, t_2, t_3)\) is composed of the three switching times:

\[\min_{y \in \mathrm{R}^3} ~J(y) \quad \text{subject to} \quad h(y) = 0 ~~\text{and}~~ y \in Y,\]

where the cost is given by

\[J(y) \overset{\mathrm{def}}{=} -\pi_B \big( \Theta(y) \big),\]

where the nonlinear equality constraints are given by

\[\begin{split}h(y) \overset{\mathrm{def}}{=} h(t_1,t_2,t_3) = \begin{pmatrix} \pi_S \big( \Psi(t_1,t_2) \big) - S^\star \\ \pi_V \big( \Theta(t_1,t_2,t_3) \big) - \bar{V} \end{pmatrix}\end{split}\]

and where the linear inequality constraints are given by

\[Y \overset{\mathrm{def}}{=} \{ y=(t_1, t_2, t_3) \in [0,1]^3 ~|~ t_1 \le t_2 \le t_3 \}.\]

The cost is simply \(-B(T)\). The first condition in \(h\) is \(S(t_2) = S^\star\) while the second one is \(V(T) = \bar{V}\). The inequality constraints are clear.

Remark. Note that since \(S(t_2) = S^\star\) and \(u(t) = u^\star(t)\) for \(t\in(t_2,t_3]\) we impose that \(S(t) = S^\star\) for every \(t\in(t_2,t_3]\).

[4]:
# Psi
def psi(t1, t2):
    x1 = exponential(Fmin, t1, t0, x0_)
    x2 = exponential(Fmax, t2, t1, x1)
    return x2

# Theta: we return also x2 for convenience
def theta(t1, t2, t3):
    x2 = psi(t1, t2)
    x3 = exponential(Fstar, t3, t2, x2)
    xf = exponential(Fmin,   T, t3, x3)
    return xf, x2

# get_states: returns B(T), V(T) and S(t2)
def get_states(t1, t2, t3):
    xf, x2 = theta(t1, t2, t3)
    return xf[1], xf[2], x2[0]

# Cost function J
def J(y):
    t1 = y[0]
    t2 = y[1]
    t3 = y[2]
    BT, VT, St2 = get_states(t1, t2, t3)
    return -BT

# Equality nonlinear constraints h
def h(y):
    t1 = y[0]
    t2 = y[1]
    t3 = y[2]
    BT, VT, St2 = get_states(t1, t2, t3)
    return np.array([VT-Vbar, St2-Ss])

nlcons = NonlinearConstraint(h, np.array([0., 0.]), np.array([0., 0.]))

# Inequality linear constraints: Y = { y in [0,1]^3 such that (0, 0) <= Ay = (t2-t1, t3-t2) }
A      = np.array([[-1., 1., 0.], [0., -1., 1.]])
lcons  = LinearConstraint(A, np.array([0., 0.]), np.array([np.inf, np.inf]))
bnds   = ((0., 1.), (0., 1.), (0., 1.))

We solve the optimization problem with scipy.optimize.minimize solver. We use the solution from [Boumaza et al. 2020] to initialize the solver.

[5]:
Nfeval = 1
print('\n Niter  t1          t2          t3          J              h')
def callbackF(yi):
    global Nfeval
    print('{0:4d}   {1: 3.6f}   {2: 3.6f}   {3: 3.6f}   {4: 3.6e}  {5: 3.6e}'.format(\
        Nfeval, yi[0], yi[1], yi[2], J(yi), np.linalg.norm(h(yi))))
    Nfeval += 1

y_init = np.array([0.6, 0.65, 0.95] ) # Initialization guessed from Figure 5 of [Boumaza et al. 2020]

res    = minimize(J, y_init, method='SLSQP', constraints={lcons, nlcons}, bounds=bnds,
                  options={'disp': True, 'ftol':1e-6, 'eps': 1e-8, 'maxiter': 100}, callback=callbackF)
ysol   = res.x

print('')
print('Optimal switching times = [{0:3.4f}, {1:3.4f}, {2:3.4f}]'.format(res.x[0], res.x[1], res.x[2]))
print('Final biomass           = {0: 3.6e}'.format(-res.fun))
print('Equality NL constraints = {0: 3.6e}'.format(np.linalg.norm(h(res.x))))

 Niter  t1          t2          t3          J              h
   1    0.535175    0.567080    0.910165   -3.605162e-02   3.544702e-03
   2    0.599855    0.634164    0.947245   -3.890431e-02   1.513621e-04
   3    0.597878    0.632127    0.945999   -3.891319e-02   2.399606e-06
   4    0.591720    0.625727    0.942570   -3.889743e-02   2.213110e-05
   5    0.596241    0.630425    0.945089   -3.891550e-02   1.621926e-08
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -0.038915495061515225
            Iterations: 5
            Function evaluations: 28
            Gradient evaluations: 5

Optimal switching times = [0.5962, 0.6304, 0.9451]
Final biomass           =  3.891550e-02
Equality NL constraints =  1.621926e-08
[6]:
# Plot the solution from direct shooting
direct = plottings.Direct(t0, x0_, T, Fmin, Fmax, Fstar, umin, umax, ustar, Sh, Sw, Ss)

t1 = ysol[0]
t2 = ysol[1]
t3 = ysol[2]
direct.plot(t1, t2, t3)
../_images/irrigation_hybrid_hybrid_15_0.png

Remark on the non differentiability of the dynamics. First, note that we use an adaptive step length Runge-Kutta scheme of high order for the numerical integration of the different systems. One can see on the following figure that around a regular crossing time (see [Haberkorn and Trélat 2011] for the definition), the numerical integrator makes shorter steps. In the following figure, we only show the first arc between \(t_0=0\) and \(t_1\approx 0.5962.\) On the figure is represented by red vertical lines, the two regular switching times, the first crossing is when \(S=S^\star\) while the second one is when \(S=S_w\). To improve efficiency (and accuracy), we should fix the system (that is the vector field) until we detect a crossing. Then, we should compute the exact crossing time and change the system accordingly to pursue the integration. See [Martinon and Gergaud 2007, Cots 2012] for details about how to implement this. We choose not to do this since we can get enough accuracy also by changing the numerical tolerances of the integrator (if needed), and because the efficiency is not an issue since there are very few crossings.

[7]:
# Plot steps length of the optimal solution for t in [t0, t1]
steps = plottings.Steps(t0, t1, Fmin, x0_, Sw, Ss)
steps.plot(Tol=1e-12)
../_images/irrigation_hybrid_hybrid_17_0.png

Main difficulty

We recall that the main difficulty is due to the fact that the optimal trajectory stays on \(M^0\) and so we have \(S(t)=S^\star\) on a time interval of non-empty interior, recalling that \(K_S\) is non differentiable at \(S^\star\). The optimal trajectory has a semi-transverse contact when reaching the stratum \(M^0\) and the associated time \(t_2\) of contact is not a regular crossing time. The optimal trajectory has another semi-transverse contact when leaving \(M^0\) at time \(t_3\) and it is neither a regular crossing time.

Back to top

Regularization and continuation

Goal 2. Retrieve, by regularization and continuation, the fact the singular arc of the optimal solution is located at a corner point of the dynamics. This is the main goal.

Regularization of the dynamics

We define the function \(K_S^{\mathrm{reg}}(S, \eta)\) such that \(K_S(S) = K_S^{\mathrm{reg}}(S, 0)\) and such that \(S \mapsto K_S^{\mathrm{reg}}(S, \eta)\) is smooth in the neighborhood of \(S^\star\) whenever \(\eta \ne 0\), \(\eta\) close to \(0\).

[8]:
# Regularization of KS
def KS_reg(S, eta):
    theta = np.arctan(Ss-Sw)
    Sc    = Sw+(Ss-Sw)*(1.0-eta+eta*np.sin(theta))+eta*np.cos(theta)
    Scm   = Sc-eta*np.cos(theta)
    if (S>Scm and S<Sc):
        KS = 1.0-eta+np.sqrt(eta**2-(S-Sc)**2)
    else:
        if(S <= Sw):
            KS = 0.
        elif(S>=Ss):
            KS = 1.
        else:
            KS = (S-Sw)/(Ss-Sw)
    return KS, Sc, Scm

# Plot the graph of KS regularized
transpi = plottings.Transpiration(Sw, Ss, KS_reg)
transpi.plot()
../_images/irrigation_hybrid_hybrid_21_0.png

The regularized problem

Let \(\eta \ne 0\), \(\eta\) close to \(0\). The regularized optimal control problem, denoted \((\mathrm{OCP}_\eta)\), is the following

\[\begin{split}\left\{ \begin{array}{l} \min J(x(\cdot), u(\cdot)) \overset{\mathrm{def}}{=} g(x(T)) = -B(T) \\[1.0em] \dot{x}(t) = F_0(t, x(t), \eta) + u(t)\, F_1, \quad u(t) \in [0, 1], \quad t \in [0, T] \text{ a.e.}, \\[1.0em] x(0) = x_0, \quad c(x(T)) = 0, \end{array} \right.\end{split}\]

where the initial state is \(x_0 \overset{\mathrm{def}}{=} (0, B_0, 0)\), where the final condition is given by \(c(x) \overset{\mathrm{def}}{=} \pi_V(x) - \bar{V}\) and where the parameterized vector field \(F_0(t, x, \eta)\) is simply the old \(F_0(t, x)\) where \(K_S\) is replaced by \(K_S^{\mathrm{reg}}(\cdot, \eta)\), that is

\[\begin{split}F_0(t, x, \eta) = \begin{pmatrix} -k_1\big(\varphi(t)\, K_S^{\mathrm{reg}}(S, \eta) + (1-\varphi(t))\, K_R(S) \big) \\[0.5em] \varphi(t)\, K_S^{\mathrm{reg}}(S, \eta)\, f(B) \\[0.5em] 0 \end{pmatrix}, \quad x=(S, B, V).\end{split}\]

Remark. The dynamics of this optimal control problem is not smooth. We only have regularized the dynamics around the corner \(S^\star\), which is sufficient for our purpose.

Application of the Pontryagin Maximum Principle

Let \((x(\cdot), u(\cdot))\) being a solution of \((\mathrm{OCP}_\eta)\) and assume that the set \(\{ t \in [0, T] ~|~ \pi_S(x(t)) = S_w \}\) is finite. Then, according to the Pontryagin Maximum Principle, there exists a covector \(p\) (absolutely continuous), a scalar \(p^0 \in \{-1, 0\}\) and a Lagrange multiplier \(\lambda \in \mathbb{R}\) such that:

  • we have the non-triviality condition:

\[(p, p^0) \ne (0, 0),\]
  • we have the pseudo-Hamiltonian differential equation

\[\dot{x}(t) = \frac{\partial H}{\partial p}(t, x(t), p(t), u(t)), \quad \dot{p}(t) = -\frac{\partial H}{\partial x}(t, x(t), p(t), u(t)),\]
  • we have the maximization condition

\[H(t, x(t), p(t), u(t)) = \max_{w \in [0, 1]} H(t, x(t), p(t), w) = H_0(t, x(t), p(t)) + \max_{w \in [0, 1]} w \, H_1(t, x(t), p(t)),\]
  • and we have the transversality condition

\[p(T) = p^0 g'(x(T)) + \lambda c'(x(T)) = (0, -p^0, \lambda),\]

where

\[H(t, x, p, u) \overset{\mathrm{def}}{=} p \cdot \left(F_0(t, x, \eta) + u\, F_1 \right) = H_0(t, x(t), p(t)) + u\, H_1(t, x(t), p(t))\]

is the pseudo-Hamiltonian associated to \((\mathrm{OCP}_\eta)\) and where \(H_0\) and \(H_1\) denote the Hamiltonian lifts of \(F_0\) and \(F_1\). Note that \(H_1\) depends only on the variable \(p\) so we could write \(H_1(p)\).

Assumption. In the following, we assume that \(p^0 = -1\), that is that we are in the normal case.

The maximization condition reads

\[\begin{split}u(t) ~ \left\{ \begin{array}{lcl} ~= -1 & \text{if} & H_1(t, x(t), p(t)) < 0, \\[0.5em] ~= +1 & \text{if} & H_1(t, x(t), p(t)) > 0, \\[0.5em] ~\in [-1, 1] & \text{if} & H_1(t, x(t), p(t)) = 0. \\[0.5em] \end{array} \right.\end{split}\]

However, if \(H_1(t, x(t), p(t)) = 0\) along a time interval on non empty interior, then we have also

\[\frac{\mathrm{d}}{\mathrm{d} t} H_1(t, x(t), p(t)) = 0.\]

Writting \(z \overset{\mathrm{def}}{=} (x,p)\), introducing the pseudo-Hamiltonian system \(\vec{H} \overset{\mathrm{def}}{=} (\partial_p H, -\partial_x H)\) and using the Poisson bracket notation, then we have

\[\begin{split}\begin{aligned} \frac{\mathrm{d}}{\mathrm{d} t} H_1(t, z(t)) &= \frac{\partial H_1}{\partial t}(t, z(t)) + \frac{\partial H_1}{\partial z}(t, z(t)) \cdot \dot{z}(t) = \frac{\partial H_1}{\partial z}(t, z(t)) \cdot \dot{z}(t) \\[0.5em] & \overset{\mathrm{def}}{=} \{H, H_1\}(t, z(t)) = \{H_0, H_1\}(t, z(t)) + u(t)\, \{H_1, H_1\}(t, z(t)) = \{H_0, H_1\}(t, z(t)) \\[0.5em] & \overset{\mathrm{def}}{=} H_{01}(t, z(t)) = 0. \end{aligned}\end{split}\]

Derivating twice, we have

\[\begin{split}\begin{aligned} \frac{\mathrm{d}}{\mathrm{d} t} H_{01}(t, z(t)) &= \frac{\partial H_{01}}{\partial t}(t, z(t)) + \frac{\partial H_{01}}{\partial z}(t, z(t)) \cdot \dot{z}(t) \\[0.5em] & = \frac{\partial H_{01}}{\partial t}(t, z(t)) + H_{001}(t, z(t)) + u(t) \, H_{101}(t, z(t)) = 0. \end{aligned}\end{split}\]

Finally, if \(H_1(t, x(t), p(t)) = 0\) along an arc of non empty interior, then the control must satisfy

\[u(t) = -\frac{\partial_t H_{01}(t, z(t)) + H_{001}(t, z(t))}{H_{101}(t, z(t))},\]

whenever \(H_{101}(t, z(t)) \ne 0\). Such an arc is called a singular arc. We define the singular control as

\[u_s(t, z) \overset{\mathrm{def}}{=} -\frac{\partial_t H_{01}(t, z) + H_{001}(t, z)}{H_{101}(t, z)}\]

such that along a singular arc, the optimal control is equal to the singular control, that is \(u(t) = u_s(t, z(t))\).

[9]:
options = nt.ivp.Options(SolverMethod='dopri5', TolAbs=1e-8, TolRel=1e-8)

# Singular control: using(t, x, p, eta)
using = wrappers.singular_control()

# Hamiltonian and flows
# The Hamiltonians are written in Fortran
# The Hamiltonian codes are compiled if docompile is True

# Hmin = H0
hmin = wrappers.hamiltonian_umin()
fmin = nt.ocp.Flow(hmin, options=options)

# Hmax = H0 + H1
hmax = wrappers.hamiltonian_umax()
fmax = nt.ocp.Flow(hmax, options=options)

# Hsing = H0 + using H1
hsing = wrappers.hamiltonian_using()
fsing = nt.ocp.Flow(hsing, options=options)

# Switching function and derivatives
h1  = wrappers.switching_function() # H1
h01 = wrappers.lift_H01()           # H01

The indirect shooting method

[10]:
# Auxiliary function
co = lambda v : np.concatenate(v)

# Shooting function: recall that he solution follows Hmin, Hmax, Hsing, Hmin
def shoot(y, eta):
    # y = (S(T), B(T), pV(T), t1, t2, t3)
    ST  = y[0]
    BT  = y[1]
    pVT = y[2]
    t1  = y[3]
    t2  = y[4]
    t3  = y[5]
    ind = 6
    nz  = 6
    z1_ = y[ind:ind+nz]; ind=ind+nz
    z2_ = y[ind:ind+nz]; ind=ind+nz
    z3_ = y[ind:ind+nz]; ind=ind+nz
    s   = np.zeros(6+3*nz)

    # xT, pT
    xT = np.zeros(3)
    pT = np.zeros(3)
    xT[0] = ST
    xT[1] = BT
    xT[2] = Vbar
    pT[0] = 0.0
    pT[1] = -pzero
    pT[2] = pVT

    x3, p3 = fmin (T,  xT, pT, t3, eta) # Hmin:  [t3, T]
    x2, p2 = fsing(t3, x3, p3, t2, eta) # Hsing: [t2, t3]
    x1, p1 = fmax (t2, x2, p2, t1, eta) # Hmax:  [t1, t2]
    x0, p0 = fmin (t1, x1, p1, t0, eta) # Hmin:  [t0, t1]

    # shooting condition
    # x(0) = x0, H1(t1)=H1(t2)=H01(t2)=0
    s[0:3] = x0 - x0_
    s[3]   = h1(p1)
    s[4]   = h1(p2)
    s[5]   = h01(t2, x2, p2, eta)
    ind = 6
    s[ind:ind+nz] = z1_ - co((x1, p1)); ind=ind+nz
    s[ind:ind+nz] = z2_ - co((x2, p2)); ind=ind+nz
    s[ind:ind+nz] = z3_ - co((x3, p3)); ind=ind+nz

    return s

# derivative of the shooting function: dshoot = shoot'(y, eta) • (dy, deta)
def dshoot(y, dy, eta, deta):
    # y = (S(T), B(T), pV(T), t1, t2, t3)
    ST, dST   = y[0], dy[0]
    BT, dBT   = y[1], dy[1]
    pVT, dpVT = y[2], dy[2]
    t1, dt1   = y[3], dy[3]
    t2, dt2   = y[4], dy[4]
    t3, dt3   = y[5], dy[5]
    ind = 6
    nz  = 6
    z1_, dz1_ = y[ind:ind+nz], dy[ind:ind+nz]; ind=ind+nz
    z2_, dz2_ = y[ind:ind+nz], dy[ind:ind+nz]; ind=ind+nz
    z3_, dz3_ = y[ind:ind+nz], dy[ind:ind+nz]; ind=ind+nz
    s, ds     = np.zeros(6+3*nz), np.zeros(6+3*nz)

    # xT, pT
    xT, dxT = np.zeros(3), np.zeros(3)
    pT, dpT = np.zeros(3), np.zeros(3)
    xT[0], dxT[0] = ST, dST
    xT[1], dxT[1] = BT, dBT
    xT[2], dxT[2] = Vbar, 0.0
    pT[0], dpT[0] = 0.0, 0.0
    pT[1], dpT[1] = -pzero, 0.0
    pT[2], dpT[2] = pVT, dpVT

    (x3, dx3), (p3, dp3) = fmin (  T,       (xT, dxT), (pT, dpT), (t3, dt3), (eta, deta)) # [t3, T]
    (x2, dx2), (p2, dp2) = fsing((t3, dt3), (x3, dx3), (p3, dp3), (t2, dt2), (eta, deta)) # [t2, t3]
    (x1, dx1), (p1, dp1) = fmax ((t2, dt2), (x2, dx2), (p2, dp2), (t1, dt1), (eta, deta)) # [t1, t2]
    (x0, dx0), (p0, dp0) = fmin ((t1, dt1), (x1, dx1), (p1, dp1),  t0,       (eta, deta)) # [t0, t1]

    # shooting condition
    # x(0) = x0, H1(t1)=H1(t2)=H01(t2)=0
    s[0:3], ds[0:3] = x0 - x0_, dx0  # x0 is a parameter
    s[3],   ds[3]   = h1((p1, dp1)) # H1 does not depend on t, x, eta
    s[4],   ds[4]   = h1((p2, dp2)) # H1 does not depend on t, x, eta
    s[5],   ds[5]   = h01((t2, dt2), (x2, dx2), (p2, dp2), (eta, deta))
    ind = 6
    s[ind:ind+nz], ds[ind:ind+nz] = z1_ - co((x1, p1)), dz1_ - co((dx1, dp1)); ind=ind+nz
    s[ind:ind+nz], ds[ind:ind+nz] = z2_ - co((x2, p2)), dz2_ - co((dx2, dp2)); ind=ind+nz
    s[ind:ind+nz], ds[ind:ind+nz] = z3_ - co((x3, p3)), dz3_ - co((dx3, dp3)); ind=ind+nz

    return s, ds

shoot = nt.tools.tensorize(dshoot, tvars=(1, 2), full=True)(shoot)
[11]:
# Shooting
if eta_init == 0.3:
    y_guess = np.array([ \
        6.05191975e-01,  3.18782899e-02, -1.00652695e+00,  5.81896245e-01,
        6.09077273e-01,  8.71129844e-01,  3.77856275e-01,  5.20074897e-04,
        1.57957791e-24,  9.58597094e-02,  5.93724551e+01, -1.00652695e+00,
        6.41091678e-01,  5.37580596e-04,  2.71810288e-02,  9.58597094e-02,
        5.74400649e+01, -1.00652695e+00,  8.30487738e-01,  3.70457096e-03,
        8.33333333e-02,  9.58597092e-02,  8.36178351e+00, -1.00652695e+00])
else:
    pass # not working

nlefun   = lambda y: shoot(y, eta_init)
sol_nle  = nt.nle.solve(nlefun, y_guess, df=nlefun)

     Calls  |f(x)|                 |x|

         1  1.555203258233392e-05  8.307548228861798e+01
         2  8.166183742856175e-06  8.307544625490236e+01
         3  9.620750526808722e-07  8.307536863186581e+01
         4  3.268830629647920e-09  8.307536826365225e+01
         5  5.882436513003873e-10  8.307536830054848e+01

 Results of the nle solver method:

 xsol    =  [ 6.05195223e-01  3.18782920e-02 -1.00649294e+00  5.81904833e-01
  6.09086627e-01  8.71130206e-01  3.77852724e-01  5.20074900e-04
  3.78858258e-24  9.58564702e-02  5.93724602e+01 -1.00649294e+00
  6.41095499e-01  5.37582332e-04  2.71817938e-02  9.58564702e-02
  5.74398848e+01 -1.00649294e+00  8.30491913e-01  3.70453448e-03
  8.33333333e-02  9.58564718e-02  8.36186622e+00 -1.00649294e+00]
 f(xsol) =  [-2.60114152e-12 -1.63139901e-15 -1.63597434e-13 -2.31547670e-10
 -1.30972344e-10 -4.24801111e-10  5.59441382e-13  1.89236647e-15
  1.63597434e-13  2.20521518e-11 -2.16040519e-10  0.00000000e+00
  5.15143483e-13  2.04111901e-15  1.63806468e-13  1.24735500e-11
 -2.17866614e-10  0.00000000e+00  5.44009282e-15  3.29597460e-17
  0.00000000e+00  4.14822343e-12  8.02913291e-13  0.00000000e+00]
 nfev    =  5
 njev    =  1
 status  =  1
 success =  True

 Successfully completed: relative error between two consecutive iterates is at most TolX.

[12]:
# Plot the solution from indirect shooting
indirect = plottings.Indirect(t0, x0_, T, umin, umax, using, fmin, fmax, fsing,
                              Sh, Sw, Ss, KS_reg, Vbar, pzero, h1, data)
indirect.plot(sol_nle.x, eta_init)
../_images/irrigation_hybrid_hybrid_28_0.png

The continuation method and final results

[13]:
# Continuation
if data.contains('path') and not restart:
    path = data.get('path');  # print('path loaded')
    eout = np.array(path.get('eout'))
    yout = np.array(path.get('yout'))
else:
    y0 = sol_nle.x
    eout, yout, nout, niter, flag = continuation.solve(shoot, y0, eta_span, df=shoot)
    if flag==1:
        print("Continuation has finished with success")
        path = {'eout':eout.tolist(),
                'yout':yout.tolist(),
                'nout':nout.tolist(),
                'niter':niter,
                'flag':flag}
        data.update({'path':path}); print('path saved')

# get final solution and print result
eta = eout[-1]
y   = yout[:,-1]
final_solution = indirect.compute_solution(y, eta)
states = np.array(final_solution['states'])

print('')
print('Solution for eta = {0:3.3f}'.format(eta))
print('')
print('  Optimal switching times = [{0:3.4f}, {1:3.4f}, {2:3.4f}]'.format(final_solution['t1'],
                                                                        final_solution['t2'],
                                                                        final_solution['t3']))
print('  Final biomass           = {0: 3.6e}'.format(states[-1,1]))
print('')

Solution for eta = 0.001

  Optimal switching times = [0.5963, 0.6304, 0.9448]
  Final biomass           =  3.889541e-02

[14]:
# Display the deformation of the solution wrt eta
indirect.interact(yout, eout, embed=True, restart=False)
[14]:
Solutions

Back to top

References

[Bock and Plitt 1984] H. G. Bock and K. J. Plitt, A Multiple Shooting Algorithm for Direct Solution of Optimal Control Problems, IFAC Proceedings Volumes, 17 (1984), no. 2, 1603-1608. Proceedings of the 9th IFAC World Congress: A Bridge Between Control Science and Technology, Budapest, Hungary, 2-6 July 1984.

[Bonnard and Chyba 2003] B. Bonnard and M. Chyba, Singular trajectories and their role in control theory, vol 40 of Mathematics & Applications, Springer-Verlag, Berlin (2003), 357 pages.

[Bonnard et al. 2015] B. Bonnard, M. Claeys, O. Cots and P. Martinon, Geometric and numerical methods in the contrast imaging problem in nuclear magnetic resonance, Acta Appl. Math., 135 (2015), no. 1, 5-45.

[Bonnard and Cots 2014] B. Bonnard and O. Cots, Geometric numerical methods and results in the control imaging problem in nuclear magnetic resonance, Math. Models Methods Appl. Sci., 24 (2014), no. 1, 187-212.

[Bonnard et al. 2020] B. Bonnard, O. Cots, J. Rouot and T. Verron, Time minimal saturation of a pair of spins and application in Magnetic Resonance Imaging, Math. Control Relat. Fields, 10 (2020), no. 1, 47-88.

[Boumaza et al. 2020] K. Boumaza, N. Kalboussi, A. Rapaport, S. Roux and C. Sinfort. Optimal control of a crop irrigation model under water scarcity, (2020) hal-01991296.

[Caines et al. 2006] P. E. Caines, F. H. Clarke, X. Liu and R. B. Vinter, A Maximum Principle for Hybrid Optimal Control Problems with Pathwise State Constraints, Proceedings of the 45th IEEE Conference on Decision and Control (2006), 4821-4825.

[Clarke 2013] F. Clarke, Functional Analysis, Calculus of Variations and Optimal Control, Springer (2013), 591 pages.

[Cots 2012] O. Cots, Contrôle optimal géométrique : méthodes homotopiques et applications. Phd thesis, Institut Mathematiques de Bourgogne, Dijon, France, 2012.

[Cots 2017] O. Cots, Geometric and numerical methods for a state constrained minimum time control problem of an electric vehicle, ESAIM Control Optim. Calc. Var., 23 (2017), no. 4, 1715-1749.

[Cots et al. 2018] O. Cots, J. Gergaud and D. Goubinat, Direct and indirect methods in optimal control with state constraints and the climbing trajectory of an aircraft, Optim. Control Appl. Meth., 39 (2018), no. 1, 281-301.

[Dmitruk and Kaganovich 2008] A. V. Dmitruk and A. M. Kaganovich, The hybrid maximum principle is a consequence of pontryagin maximum principle, Systems & Control Letters, 57 (2008), no. 11, 964-970.

[Haberkorn and Trélat 2011] T. Haberkorn and E. Trélat, Convergence results for smooth regularizations of hybrid nonlinear optimal control problems, SIAM J. Control Optim., 49 (2011), no. 4, 1498-1522.

[Martinon and Gergaud 2007] P. Martinon and J. Gergaud, Using switching detection and variational equations for the shooting method, Optimal Control Appl. Methods, 28 (2007), no. 2, 95–116.

[Sussmann 1999] H. J. Sussmann, A nonsmooth hybrid maximum principle, in Stability and Stabilization of Nonlinear Systems (Ghent, 1999), Lecture Notes in Control and Inform. Sci. 246, Springer, London (1999), 325–354.

[Sussmann 1999 (GOC)] H. J. Sussmann, Geometry and Optimal Control, in Baillieul J., Willems J.C. (eds) Mathematical Control Theory, Springer, New York, NY (1999), 140-198.

Back to top

Annexes

Hamiltonians lifts: h0fun.f90, h1fun.f90

[15]:
!pygmentize h0fun.f90
print('')
!pygmentize h1fun.f90
Subroutine h0fun(t,n,x,p,npar,par,h0)
    implicit none
    integer, intent(in)                             :: n, npar
    double precision, intent(in)                    :: t
    double precision, dimension(n),    intent(in)   :: x, p
    double precision, dimension(npar), intent(in)   :: par
    double precision, intent(out)                   :: h0

    !Local declarations
    double precision    :: KS, KR, phi, f
    double precision    :: S, B, theta, Sc, Scmoins
    double precision    :: Tf, k1, k2, Ss, Sw, Sh, Fmax, Qbar, alpha, r, Bmax, B0, eta, pzero

    ! Parameters
    call parameters(npar, par, Tf, k1, k2, Ss, Sw, Sh, Fmax, Qbar, alpha, r, Bmax, B0, eta, pzero)

    ! Variables
    S = x(1)
    B = x(2)
    
    ! KS Regularization
    theta   = atan(Ss-Sw)
    Sc      = Sw+(Ss-Sw)*(1d0-eta+eta*sin(theta))+eta*cos(theta)
    Scmoins = Sc-eta*cos(theta)
    if ((S>Scmoins) .and. (S<Sc)) then
        KS = 1d0-eta+sqrt(eta**2-(S-Sc)**2)
    else
        if(S.le.Sw)then
            KS = 0d0
        elseif(S.ge.Ss)then
            KS = 1d0
        else
            KS = (S-Sw)/(Ss-Sw)
        end if
    endif

    ! KR
    if(S.le.Sh)then
        KR = 0d0
    elseif(S.ge.1d0)then
        KR = 1d0
    else
        KR = (S-Sh)/(1d0-Sh)
    end if

    ! phi and f
    phi = (t/Tf)**alpha
!    f   = 1d0 !r*B*(1d0-B/Bmax)
    f   = r*B*(1d0-B/Bmax)
    
    ! H0
    h0 = -p(1)*k1*(phi*KS+(1d0-phi)*KR) + p(2)*phi*KS*f

end subroutine h0fun

Subroutine h1fun(n,p,npar,par,h1)
    implicit none
    integer, intent(in)                             :: n,npar
    double precision, dimension(n),    intent(in)   :: p
    double precision, dimension(npar), intent(in)   :: par
    double precision, intent(out)                   :: h1

    !Local declarations
    double precision    :: Tf, k1, k2, Ss, Sw, Sh, Fmax, Qbar, alpha, r, Bmax, B0, eta, pzero

    ! Parameters
    call parameters(npar, par, Tf, k1, k2, Ss, Sw, Sh, Fmax, Qbar, alpha, r, Bmax, B0, eta, pzero)

    ! H1
    h1 = p(1)*k1*k2 + p(3)

end subroutine h1fun

plottings.py

[16]:
!pygmentize plottings.py
import numpy as np
import nutopy as nt

import os.path
import panel as pn
css = '''
.bk.panel-widget-box {
  background: #f0f0f0;
  border-radius: 5px;
  border: 1px black solid;
}
'''
pn.extension('katex', 'mathjax', raw_css=[css])

from IPython.display import display, HTML
from matplotlib.figure import Figure

import matplotlib.pyplot as plt

from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes
from mpl_toolkits.axes_grid1.inset_locator import mark_inset

# colors
blue  = '#1f77b4'
green = '#2ca02c'
red   = '#d62728'

class Indirect:

    def __init__(self, t0, x0, T, umin, umax, using, fmin, fmax, fsing, Sh, Sw, Ss, KS_regul, Vbar, pzero, h1, data):
        # parameters
        # t0, x0, T, fmin, fmax, fsing, umin, umax, using, Sh, Sw, Ss, KS_regul
        self.__t0    = t0
        self.__x0    = x0
        self.__T     = T
        self.__fmin  = fmin
        self.__fmax  = fmax
        self.__fsing = fsing
        self.__umin  = umin
        self.__umax  = umax
        self.__using = using
        self.__Sh    = Sh
        self.__Sw    = Sw
        self.__Ss    = Ss
        self.__KS_regul = KS_regul
        self.__Vbar = Vbar
        self.__pzero = pzero
        self.__h1 = h1
        #
        # other
        self.__solutions_key = 'indirect_solutions'
        self.__data = data
        self.fig        = None
        self.axes       = None
        self.embeded    = False

    def __get_parameters__(self):
        # parameters
        t0    = self.__t0
        x0    = self.__x0
        T     = self.__T
        fmin  = self.__fmin
        fmax  = self.__fmax
        fsing = self.__fsing
        umin  = self.__umin
        umax  = self.__umax
        using = self.__using
        Sh    = self.__Sh
        Sw    = self.__Sw
        Ss    = self.__Ss
        KS_regul = self.__KS_regul
        Vbar = self.__Vbar
        pzero = self.__pzero
        h1 = self.__h1
        return t0, x0, T, fmin, fmax, fsing, umin, umax, using, Sh, Sw, Ss, KS_regul, Vbar, pzero, h1

    def __decorate__(self, ax, num, interactive):

        # parameters
        t0, x0, T, fmin, fmax, fsing, umin, umax, using, Sh, Sw, Ss, KS_regul, Vbar, pzero, h1 = self.__get_parameters__()

        font_size = 12

        dx = 0.05
        gap = 5e-2

        # B
        if num==0:
            ax.set_xlabel('t');
            ax.set_title('Biomass (B)')
            ax.axvline(T,  color='k', linestyle='--' , linewidth=1.);
            ax.spines['right'].set_visible(False); ax.spines['top'].set_visible(False) # remove some axes
            ax.plot((1), (0), ls="", marker=">", ms=5, color="k", transform=ax.get_yaxis_transform(), clip_on=False) # make arrows
            ax.plot((0), (1), ls="", marker="^", ms=5, color="k", transform=ax.get_xaxis_transform(), clip_on=False)
            if interactive:
                for item in ([ax.title, ax.xaxis.label, ax.yaxis.label] + ax.get_xticklabels() + ax.get_yticklabels()):
                    item.set_fontsize(font_size)
            ax.set_xlim(0.0, 1.0+gap)
            ax.set_ylim(0.0, 0.04+1e-3)

        # u
        if num==1:
            ax.set_xlabel('t');
            ax.set_title('Optimal control (u)')
            ax.axhline(1,  color='k', linestyle='--' , linewidth=1.);
            ax.axvline(T,  color='k', linestyle='--' , linewidth=1.);
            ax.spines['right'].set_visible(False); ax.spines['top'].set_visible(False) # remove some axes
            ax.plot((1), (0), ls="", marker=">", ms=5, color="k", transform=ax.get_yaxis_transform(), clip_on=False) # make arrows
            ax.plot((0), (1), ls="", marker="^", ms=5, color="k", transform=ax.get_xaxis_transform(), clip_on=False)
            if interactive:
                for item in ([ax.title, ax.xaxis.label, ax.yaxis.label] + ax.get_xticklabels() + ax.get_yticklabels()):
                    item.set_fontsize(font_size)
            ax.set_xlim(0.0, 1.0+gap)
            ax.set_ylim(0.0, 1.0+gap)

        # S
        if num==2:
            ax.set_xlabel('t');
            ax.set_title('Humidity (S)')
            ax.axhline(1,  color='k', linestyle='--' , linewidth=1.);
            ax.axvline(T,  color='k', linestyle='--' , linewidth=1.);
            ax.axhline(Ss, color='k', linewidth=1., linestyle='--' );
            ax.axhline(Sw, color='k', linewidth=1., linestyle='--' );
            ax.axhline(Sh, color='k', linewidth=1., linestyle='--' );
            ax.spines['right'].set_visible(False); ax.spines['top'].set_visible(False) # remove some axes
            ax.plot((1), (0), ls="", marker=">", ms=5, color="k", transform=ax.get_yaxis_transform(), clip_on=False) # make arrows
            ax.plot((0), (1), ls="", marker="^", ms=5, color="k", transform=ax.get_xaxis_transform(), clip_on=False)
            if interactive:
                for item in ([ax.title, ax.xaxis.label, ax.yaxis.label] + ax.get_xticklabels() + ax.get_yticklabels()):
                    item.set_fontsize(font_size)
            ax.set_xlim(0.0, 1.0+gap)
            ax.set_ylim(0.0, 1.0+gap)
            ax.set_yticks((0., Sh, Sw, Ss, 1.))
            ax.set_yticklabels(['0', '$S_h$', '$S_w$', '$S^\star$', '1'])

        # H1
        if num==3:
            yo = -1.4
            ax.set_xlabel('t');
            ax.set_title('Switching function $H_1$')
            ax.axhline(0,  color='k', linestyle='--' , linewidth=1.);
            ax.axvline(T,  color='k', linestyle='--' , linewidth=1.);
            ax.spines['right'].set_visible(False); ax.spines['top'].set_visible(False) # remove some axes
            ax.plot((1), (yo), ls="", marker=">", ms=5, color="k", transform=ax.get_yaxis_transform(), clip_on=False) # make arrows
            ax.plot((0), (1), ls="", marker="^", ms=5, color="k", transform=ax.get_xaxis_transform(), clip_on=False)
            if interactive:
                for item in ([ax.title, ax.xaxis.label, ax.yaxis.label] + ax.get_xticklabels() + ax.get_yticklabels()):
                    item.set_fontsize(font_size)
            ax.set_xlim(0.0, 1.0+gap)
            ax.set_ylim(yo, 0.0+gap)

        # KS
        if num==4:
            o = 0.5
            ax.set_xlabel('S');
            ax.set_title('$K_S$ and $K_S^{\mathrm{reg}}(\cdot, \eta)$')
            ax.plot([0., 1.], [1., 1.],  color='k', linestyle='--' , linewidth=1.);
            ax.plot([1., 1.], [0., 1.],  color='k', linestyle='--' , linewidth=1.);
            ax.plot([Ss, Ss], [0., 1.],  color='k', linestyle='--' , linewidth=1.);
            ax.spines['right'].set_visible(False); ax.spines['top'].set_visible(False) # remove some axes
            ax.plot((1), (o), ls="", marker=">", ms=5, color="k", transform=ax.get_yaxis_transform(), clip_on=False) # make arrows
            ax.plot((o), (1), ls="", marker="^", ms=5, color="k", transform=ax.get_xaxis_transform(), clip_on=False)
            if interactive:
                for item in ([ax.title, ax.xaxis.label, ax.yaxis.label] + ax.get_xticklabels() + ax.get_yticklabels()):
                    item.set_fontsize(font_size)
            ax.set_xticks((0., Sw, Ss, 1.))
            ax.set_xticklabels(['0', '$S_w$', '$S^\star$', '1'])
            ax.set_yticks((0., 1.))
            ax.set_yticklabels(['0', '1'])
            ax.set_xlim([o,1+gap])
            ax.set_ylim([o,1+gap]);

    def __init_plot__(self, interactive):

        # figure
        if interactive:
            fig = Figure(dpi=150)
            fig.set_figwidth(9)
            fig.set_figheight(8)
            fig.subplots_adjust(left=0.1, right=0.9, bottom=0.1, top=0.95, wspace=0.3, hspace=0.3)
        else:
            fig = plt.figure(dpi=100)
            fig.set_figwidth(9)
            fig.set_figheight(8)
            fig.subplots_adjust(left=0, right=1, bottom=0, top=1, wspace=0.3, hspace=0.3)

        # subplots
        num = 0
        ax = fig.add_subplot(231); self.__decorate__(ax, num, interactive); num = num + 1;
        ax = fig.add_subplot(232); self.__decorate__(ax, num, interactive); num = num + 1;
        ax = fig.add_subplot(233); self.__decorate__(ax, num, interactive); num = num + 1;
        ax = fig.add_subplot(223); self.__decorate__(ax, num, interactive); num = num + 1;
        ax = fig.add_subplot(224); self.__decorate__(ax, num, interactive); num = num + 1;

        return fig

    def __key__(self, eta):
        return str(hash(round(eta, 10)))

    def __get_solution__(self, y, eta):

        # if saved, then load, else compute

        # check if bdd has the dict containing the solutions
        # if not then create it and compute the solution, update the dict and add it to the bdd
        # if yes then
        #    if the sol has already been computed then get it from the bdd and return it
        #    else compute it, add it to the bdd and return it
        if self.__data.contains(self.__solutions_key):
            solutions = self.__data.get(self.__solutions_key)
            if self.__key__(eta) in solutions:
                sol = solutions[self.__key__(eta)]
            else:
                sol = self.__compute_solution__(y, eta)
                solutions[self.__key__(eta)] = sol
                self.__data.update({self.__solutions_key:solutions})
        else:
            solutions = {}
            sol = self.__compute_solution__(y, eta)
            solutions[self.__key__(eta)] = sol
            self.__data.update({self.__solutions_key:solutions})

        return self.__expand_solution__(sol)

    def __expand_solution__(self, sol):
        return np.array(sol['times']), np.array(sol['states']), np.array(sol['costates']), \
                np.array(sol['control']), np.array(sol['switching']), sol['t1'], sol['t2'], sol['t3'], \
                np.array(sol['S_sing']), np.array(sol['K_sing']), np.array(sol['Si']), np.array(sol['KSi_0']), np.array(sol['KSi'])

    def __agreg_solution__(self, times, states, costates, control, switching, t1, t2, t3, S_sing, K_sing, Si, KSi_0, KSi):
        sol = {}
        sol['times'] = times.tolist()
        sol['states'] = states.tolist()
        sol['costates'] = costates.tolist()
        sol['control'] = control.tolist()
        sol['switching'] = switching.tolist()
        sol['t1'] = t1
        sol['t2'] = t2
        sol['t3'] = t3
        sol['S_sing'] = S_sing.tolist()
        sol['K_sing'] = K_sing.tolist()
        sol['Si'] = Si.tolist()
        sol['KSi_0'] = KSi_0.tolist()
        sol['KSi'] = KSi.tolist()
        return sol

    def compute_solution(self, y, eta):
        return self.__compute_solution__(y, eta)

    def __compute_solution__(self, y, eta):

        # parameters
        t0, x0, T, fmin, fmax, fsing, umin, umax, using, Sh, Sw, Ss, KS_regul, Vbar, pzero, h1 = self.__get_parameters__()

        # get p0
        ST  = y[0]
        BT  = y[1]
        pVT = y[2]
        t1  = y[3]
        t2  = y[4]
        t3  = y[5]

        xT = np.zeros(3)
        pT = np.zeros(3)
        xT[0] = ST
        xT[1] = BT
        xT[2] = Vbar
        pT[0] = 0.0
        pT[1] = -pzero
        pT[2] = pVT

        x3, p3 = fmin (T,  xT, pT, t3, eta)   # Hmin: [t3, T]
        x2, p2 = fsing(t3, x3, p3, t2, eta) # Hsing: [t2, t3]
        x1, p1 = fmax (t2, x2, p2, t1, eta)  # Hmax: [t1, t2]
        x0, p0 = fmin (t1, x1, p1, t0, eta) # Hmin: [t0, t1]

        # compute the trajectory to plot
        N      = 50

        tspan1  = list(np.linspace(t0, t1, N+1))
        tspan2  = list(np.linspace(t1, t2, N+1))
        tspan3  = list(np.linspace(t2, t3, N+1))
        tspanf  = list(np.linspace(t3,  T, N+1))

        x1, p1 = fmin (t0, x0, p0, tspan1, eta)
        x2, p2 = fmax (t1, x1[-1], p1[-1], tspan2, eta)
        x3, p3 = fsing(t2, x2[-1], p2[-1], tspan3, eta)
        xf, pf = fmin (t3, x3[-1], p3[-1], tspanf, eta)

        u1     = umin(tspan1)
        u2     = umax(tspan2)
        u3     = using(tspan3, x3, p3, eta)
        uf     = umin(tspanf)

        times    = np.concatenate((tspan1, tspan2, tspan3, tspanf))
        states   = np.array(np.concatenate((x1, x2, x3, xf)))
        costates = np.array(np.concatenate((p1, p2, p3, pf)))
        control  = np.array(np.concatenate((u1, u2, u3, uf)))

        h11      = h1(p1) #h1(tspan1, x1, p1, eta)
        h12      = h1(p2) #h1(tspan2, x2, p2, eta)
        h13      = h1(p3) #h1(tspan3, x3, p3, eta)
        h1f      = h1(pf) #h1(tspanf, xf, pf, eta)

        switching = np.array(np.concatenate((h11, h12, h13, h1f)))

        S_sing = np.array(x3)[:,0]
        N_sing = len(S_sing)
        K_sing = np.array([(KS_regul(S_sing[i], eta)[0]) for i in range(0, N_sing)])

        Si  = np.concatenate((np.linspace(0,Sw,50), np.linspace(Sw,Ss,50), np.linspace(Ss,1.,50)))
        KSi_0 = np.zeros(Si.size)
        KSi   = np.zeros(Si.size)
        for i in range(0,Si.size):
            KSi_0[i], _ , _ = KS_regul(Si[i], 0.0)
            KSi[i], Sc , Scm   = KS_regul(Si[i], eta)

        return self.__agreg_solution__(times, states, costates, control, switching, t1, t2, t3, S_sing, K_sing, Si, KSi_0, KSi)

    # static plot
    def plot(self, y, eta, fig=None, newfig=True, interactive=False, zoom=True):
        return self.__plot__(y, eta, fig=fig, newfig=newfig, interactive=interactive, zoom=zoom)

    def __plot__(self, y, eta, fig, newfig, interactive, zoom): # function to plot a control-state-costate trajectory with structure 0, 1, us, 1

        # if fig is None, then we initialize the plot
        if interactive:
            self.fig = self.__init_plot__(interactive)
        else:
            if fig is None:
                if self.fig is None or newfig:
                    self.fig = self.__init_plot__(interactive)
            else:
                self.fig = fig

        #
        self.axes = self.fig.get_axes()
        if interactive:
            for num in range(0, len(self.axes)):
                self.axes[num].clear()
                self.__decorate__(self.axes[num], num, interactive)

        # parameters
        t0, x0, T, fmin, fmax, fsing, umin, umax, using, Sh, Sw, Ss, KS_regul, Vbar, pzero, h1 = self.__get_parameters__()

        # solution
        times, states, costates, control, switching, t1, t2, t3, S_sing, K_sing, Si, KSi_0, KSi = self.__get_solution__(y, eta)

        # color for singular arc part
        c_sing = (0.8, 0.8, 0.8, 0.5)

        # -------------------------------------------------------
        # B
        ax = self.axes[0]
        ax.plot(times, states[:,1])

        # -------------------------------------------------------
        # u
        ax = self.axes[1]
        ax.fill_between([t2, t3], [0, 0], [1, 1], color=c_sing, edgecolor=None, linewidth=0.0)
        ax.plot(times, control)

        # -------------------------------------------------------
        # S
        ax = self.axes[2]
        S_sing_min = min(S_sing)
        S_sing_max = max(S_sing)
        ax.fill_between([0, 1], [S_sing_min, S_sing_min], [S_sing_max, S_sing_max], color=c_sing, edgecolor=None, linewidth=0.0)
        ax.fill_between([t2, t3], [0, 0], [S_sing_min, S_sing_min], color=c_sing, edgecolor=None, linewidth=0.0)
        ax.fill_between([t2, t3], [S_sing_max, S_sing_max], [1, 1], color=c_sing, edgecolor=None, linewidth=0.0)
        ax.plot(times, states[:,0])

        # -------------------------------------------------------
        # H1
        ax = self.axes[3]
        yl = ax.get_ylim()
        ax.fill_between([t2, t3], [yl[0], yl[0]], [yl[1], yl[1]], color=c_sing, edgecolor=None, linewidth=0.0)
        ax.axvline(t1, color='k', linestyle='--', linewidth=0.5);
        ax.axvline(t2, color='k', linestyle='--', linewidth=0.5);
        ax.axvline(t3, color='k', linestyle='--', linewidth=0.5);
        ax.plot(times, switching)

        if zoom:
            # zoom
            axins = zoomed_inset_axes(ax, 13, loc='lower center', borderpad=1)
            axins.fill_between([t2, t3], [yl[0], yl[0]], [yl[1], yl[1]], color=c_sing, edgecolor=None, linewidth=0.0)
            axins.axvline(t1, color='k', linestyle='--', linewidth=0.5);
            axins.axvline(t2, color='k', linestyle='--', linewidth=0.5);
            axins.axvline(t3, color='k', linestyle='--', linewidth=0.5);
            axins.axhline(0., color='k', linewidth=0.5, linestyle='--' );
            axins.plot(times, switching)

            # sub region of the original image
            axins.set_xlim(0.58, 0.64)
            axins.set_ylim(-0.015, 0.015)

            axins.get_xaxis().set_visible(False)
            axins.get_yaxis().set_visible(False)

            # draw a bbox of the region of the inset axes in the parent axes and
            # connecting lines between the bbox and the inset axes area
            mark_inset(ax, axins, loc1=2, loc2=1, fc="none", ec="0.5")

        # -------------------------------------------------------
        # KS
        ax = self.axes[4]
        #ax.fill_between(S_sing, K_sing, color=c_sing, edgecolor=None)
        ax.fill_between([S_sing_min, S_sing_max], [1, 1], color=c_sing)
        ax.plot(Si, KSi_0, color=red);
        ax.plot(Si, KSi);

        if interactive:
            return self.fig

    # Interactive plot
    def interact(self, ys, etas, embed=False, restart=False):
        # compute solutions to save it if needed
        for i in range(0,len(etas)):
            y = ys[:,i]
            eta = etas[i]
            times, states, costates, control, switching, t1, t2, t3, S_sing, K_sing, Si, KSi_0, KSi = self.__get_solution__(y, eta)

        #
        filename = 'indirect_solutions.html'
        if not( embed and os.path.isfile(filename) ) or restart:
            pane = self.__plot_interactive__(ys, etas, embed=embed, filename=filename)
        else:
            pane = None
        if embed:
            return HTML(filename)
        else:
            return pane

    #
    def __plot_interactive__(self, ys, etas, embed, filename):

        def myplot(num=0):
            return self.__plot__(ys[:, num], etas[num], fig=None, newfig=False, interactive=True, zoom=True)

        # Slider haut
        num_sol = pn.widgets.Player(name='Player', start=0, end=len(etas)-1, value=0, interval=343, loop_policy='reflect',
                            margin=(0, 20, 0, 20), width_policy='max', sizing_mode='scale_both') # H, D, B, G)

        value_num_sol_1 = pn.pane.LaTeX(r"""$\eta =$""", renderer='mathjax', style={'font-size': '12pt'}, margin=(0, 0, 0, 0))
        if embed:
            m = (0, 0, 0, -25)
        else:
            m = (0, 0, 0, 5)
        value_num_sol_2 = pn.pane.LaTeX(renderer='mathjax', style={'font-size': '12pt'}, margin=m)
        value_num_sol = pn.GridBox(value_num_sol_1, value_num_sol_2, sizing_mode='scale_both', ncols=2, align='center', margin=(10, 0, 0, 0))

        # Controlleur haut
        @pn.depends(num_sol.param.value)
        def callback_num_sol(num):
            value_num_sol_2.object='{0:3.3f}'.format(etas[num])

        # Figures
        reactive_plot = pn.bind(myplot, num_sol)
        figures = pn.panel(reactive_plot, sizing_mode='scale_both', margin=(-30, 0, 0, 0), align='center')

        # Legend
        legend = pn.pane.LaTeX(r"""Figures: legend to be written.""", renderer='mathjax', style={'font-size': '12pt'}, width_policy='max', align='center')

        # Pane
        pane = pn.GridBox(callback_num_sol, figures, legend, value_num_sol, num_sol, \
                      max_width=800, sizing_mode='scale_both', align='center')

        if embed:
            if not self.embeded:
                step = 1
                pane.save(filename, embed_states={num_sol: list(range(0, len(etas), step))}, progress=True, embed=True, title='Solutions')
                self.embeded = True

        return pane

class Transpiration:

    def __init__(self, Sw, Ss, KS_regul):
        # Sw, Ss, KS_regul
        self.__Sw       = Sw
        self.__Ss       = Ss
        self.__KS_regul = KS_regul

    def plot(self):

        # colors
        blue  = '#1f77b4'
        green = '#2ca02c'
        red   = '#d62728'

        # parameters
        Sw       = self.__Sw
        Ss       = self.__Ss
        KS_regul = self.__KS_regul

        # Graphs of KS regularized
        plt.rcParams.update({'font.size':12})
        fig = plt.figure(dpi=120); fig.set_figwidth(9)
        ax  = fig.add_subplot(111);
        ax.plot([0., 1.], [1., 1.],  color='k', linestyle='--' , linewidth=1.);
        ax.plot([1., 1.], [0., 1.],  color='k', linestyle='--' , linewidth=1.);
        ax.plot([Ss, Ss], [0., 1.],  color='k', linestyle='--' , linewidth=1.);
        plt.xticks((0., Sw, Ss, 1.), ['0', '$S_w$', '$S^\star$', '1'])
        plt.yticks((0., 1.), ['0', '1'])
        plt.xlim([0,1+5e-2]); plt.ylim([0,1+5e-2]);
        ax.spines['right'].set_visible(False); ax.spines['top'].set_visible(False) # remove some axes
        ax.plot((1), (0), ls="", marker=">", ms=5, color="k", transform=ax.get_yaxis_transform(), clip_on=False) # make arrows
        ax.plot((0), (1), ls="", marker="^", ms=5, color="k", transform=ax.get_xaxis_transform(), clip_on=False)

        #
        Si  = np.concatenate((np.linspace(0,Sw,50), np.linspace(Sw,Ss,50), np.linspace(Ss,1.,50)))

        #
        eta = 0.
        KSi = np.zeros(Si.size)
        for i in range(0,Si.size):
            KSi[i], _, _ = KS_regul(Si[i], eta)
        ax.plot(Si, KSi, color=blue, linewidth=2.0, zorder=10, label='$K_S(S)$');

        #
        eta = 0.15
        KSi = np.zeros(Si.size)
        for i in range(0,Si.size):
            KSi[i], _, _ = KS_regul(Si[i], eta)
        ax.plot(Si, KSi, color=red, linewidth=2.0, zorder=10, label='$K_S^{\mathrm{reg}}(S, 0.15)$');

        #
        eta = 0.3
        KSi = np.zeros(Si.size)
        for i in range(0,Si.size):
            KSi[i], _ , _ = KS_regul(Si[i], eta)
        ax.plot(Si, KSi, color=green, linewidth=2.0, zorder=10, label='$K_S^{\mathrm{reg}}(S, 0.3)$');

        #
        ax.legend(fontsize=10, loc='upper left', bbox_to_anchor=(0.03, 0.9));


class Steps:

    def __init__(self, t0, t1, Fmin, x0, Sw, Ss):
        # t0, t1, Fmin, x0, Sw, Ss
        self.__t0    = t0
        self.__x0    = x0
        self.__t1    = t1
        self.__Fmin  = Fmin
        self.__Sw    = Sw
        self.__Ss    = Ss

    def plot(self, Tol):

        # colors
        blue  = '#1f77b4'
        green = '#2ca02c'
        red   = '#d62728'

        # parameters
        t0    = self.__t0
        x0    = self.__x0
        t1    = self.__t1
        Fmin  = self.__Fmin
        Sw    = self.__Sw
        Ss    = self.__Ss

        # Compute x for t in [t0, t1]
        opt = nt.ivp.Options(TolAbs=Tol, TolRel=Tol)
        sol = nt.ivp.exp(Fmin, t1, t0, x0, options=opt); xs = sol.xout; ts = sol.tout;

        # Get times when S = Ss and S = Sw
        istar = np.argmin(np.abs(xs[:,0]-Ss))
        iw    = np.argmin(np.abs(xs[:,0]-Sw))

        # plot step length with respect to time
        fig = plt.figure(dpi=100)
        fig.set_figwidth(9)
        ax  = fig.add_subplot(111);
        ax.axvline(t0, color='k', linewidth=0.5);
        ax.axvline(t1, color='k', linewidth=0.5);
        ax.axvline(ts[istar], color=red, linewidth=1.);
        ax.axvline(ts[iw], color=red, linewidth=1.);
        ax.set_xlabel('$t \in [t_0,t_1]$'); ax.set_ylabel('step length', rotation=0); ax.yaxis.set_label_coords(-0.1,0.5)
        plt.text(0.12, 7*1e-4, '$S=S^\star$', color=red)
        plt.text(0.46, 7*1e-4, '$S=S_w$', color=red)
        ax.plot(ts[0:-1], ts[1:]-ts[0:-1], linewidth=2.)
        plt.yscale('log')

class Direct:

    def __init__(self, t0, x0, T, Fmin, Fmax, Fstar, umin, umax, ustar, Sh, Sw, Ss):
        # t0, x0, T, Fmin, Fmax, Fstar, umin, umax, ustar, Sh, Sw, Ss
        self.__t0    = t0
        self.__x0    = x0
        self.__T     = T
        self.__Fmin  = Fmin
        self.__Fmax  = Fmax
        self.__Fstar = Fstar
        self.__umin  = umin
        self.__umax  = umax
        self.__ustar = ustar
        self.__Sh    = Sh
        self.__Sw    = Sw
        self.__Ss    = Ss

    def plot(self, t1, t2, t3):
        # function to plot a control-state trajectory with structure 0, 1, us, 0
        # t1, t2 and t3 being the three switching times

        # colors
        blue  = '#1f77b4'
        green = '#2ca02c'
        red   = '#d62728'

        # parameters
        t0    = self.__t0
        x0    = self.__x0
        T     = self.__T
        Fmin  = self.__Fmin
        Fmax  = self.__Fmax
        Fstar = self.__Fstar
        umin  = self.__umin
        umax  = self.__umax
        ustar = self.__ustar
        Sh    = self.__Sh
        Sw    = self.__Sw
        Ss    = self.__Ss

        N      = 20

        tspan1  = list(np.linspace(t0, t1, N+1))
        tspan2  = list(np.linspace(t1, t2, N+1))
        tspan3  = list(np.linspace(t2, t3, N+1))
        tspanf  = list(np.linspace(t3,  T, N+1))

        sol    = nt.ivp.exp(Fmin,  t1, t0, x0,     time_steps=tspan1); x1 = sol.xout
        sol    = nt.ivp.exp(Fmax,  t2, t1, x1[-1], time_steps=tspan2); x2 = sol.xout
        sol    = nt.ivp.exp(Fstar, t3, t2, x2[-1], time_steps=tspan3); x3 = sol.xout
        sol    = nt.ivp.exp(Fmin,   T, t3, x3[-1], time_steps=tspanf); xf = sol.xout

        u1     = umin(tspan1)
        u2     = umax(tspan2)
        u3     = ustar(tspan3)
        uf     = umin(tspanf)

        times   = np.concatenate((tspan1, tspan2, tspan3, tspanf))
        states  = np.array(np.concatenate((x1, x2, x3, xf)))
        control = np.array(np.concatenate((u1, u2, u3, uf)))

        fig = plt.figure(dpi=100)
        fig.set_figwidth(9)
        fig.set_figheight(4)
        fig.subplots_adjust(left=0, right=1, bottom=0, top=1, wspace=0.3, hspace=0.3)

        #
        gap = 5e-2
        # color for singular arc part
        c_sing = (0.8, 0.8, 0.8, 0.5)

        # B
        ax  = fig.add_subplot(131);
        ax.set_xlabel('t');
        ax.set_title('Biomass (B)')
        ax.axvline(T,  color='k', linestyle='--' , linewidth=1.);
        ax.spines['right'].set_visible(False); ax.spines['top'].set_visible(False) # remove some axes
        ax.plot((1), (0), ls="", marker=">", ms=5, color="k", transform=ax.get_yaxis_transform(), clip_on=False) # make arrows
        ax.plot((0), (1), ls="", marker="^", ms=5, color="k", transform=ax.get_xaxis_transform(), clip_on=False)
        ax.set_xlim(0.0, 1.0+gap)
        ax.set_ylim(0.0, 0.04+1e-3)
        ax.plot(times, states[:,1])

        # u
        ax  = fig.add_subplot(132);
        ax.set_xlabel('t');
        ax.set_title('Optimal control (u)')
        ax.axhline(1,  color='k', linestyle='--' , linewidth=1.);
        ax.axvline(T,  color='k', linestyle='--' , linewidth=1.);
        ax.spines['right'].set_visible(False); ax.spines['top'].set_visible(False) # remove some axes
        ax.plot((1), (0), ls="", marker=">", ms=5, color="k", transform=ax.get_yaxis_transform(), clip_on=False) # make arrows
        ax.plot((0), (1), ls="", marker="^", ms=5, color="k", transform=ax.get_xaxis_transform(), clip_on=False)
        ax.set_xlim(0.0, 1.0+gap)
        ax.set_ylim(0.0, 1.0+gap)
        ax.fill_between([t2, t3], [0, 0], [1, 1], color=c_sing, edgecolor=None, linewidth=0.0)
        ax.plot(times, control)

        ax  = fig.add_subplot(133);
        ax.set_xlabel('t');
        ax.set_title('Humidity (S)')
        ax.axhline(1,  color='k', linestyle='--' , linewidth=1.);
        ax.axvline(T,  color='k', linestyle='--' , linewidth=1.);
        ax.axhline(Ss, color='k', linewidth=1., linestyle='--' );
        ax.axhline(Sw, color='k', linewidth=1., linestyle='--' );
        ax.axhline(Sh, color='k', linewidth=1., linestyle='--' );
        ax.spines['right'].set_visible(False); ax.spines['top'].set_visible(False) # remove some axes
        ax.plot((1), (0), ls="", marker=">", ms=5, color="k", transform=ax.get_yaxis_transform(), clip_on=False) # make arrows
        ax.plot((0), (1), ls="", marker="^", ms=5, color="k", transform=ax.get_xaxis_transform(), clip_on=False)
        ax.set_xlim(0.0, 1.0+gap)
        ax.set_ylim(0.0, 1.0+gap)
        ax.set_yticks((0., Sh, Sw, Ss, 1.))
        ax.set_yticklabels(['0', '$S_h$', '$S_w$', '$S^\star$', '1'])
        ax.fill_between([t2, t3], [0, 0], [1, 1], color=c_sing, edgecolor=None, linewidth=0.0)
        ax.plot(times, states[:,0])

class KSKR:

    def __init__(self, Sh, Sw, Ss):
        # Sh, Sw, Ss
        self.__Sh = Sh
        self.__Sw = Sw
        self.__Ss = Ss

    def plot(self):

        # Colors
        blue  = '#1f77b4'
        green = '#2ca02c'
        red   = '#d62728'

        # Params
        Sh = self.__Sh
        Sw = self.__Sw
        Ss = self.__Ss

        # Init figure for graphs of KS and KR
        plt.rcParams.update({'font.size':12})
        fig = plt.figure(dpi=120);
        fig.set_figheight(3); fig.set_figwidth(7);

        # Graph of KS
        ax  = fig.add_subplot(121);
        plt.text(0.33, 0.5, '$K_S(S)$', color=blue)
        ax.plot([0., 1.], [1., 1.],  color='k', linestyle='--' , linewidth=1.);
        ax.plot([1., 1.], [0., 1.],  color='k', linestyle='--' , linewidth=1.);
        ax.plot([Ss, Ss], [0., 1.],  color=red, linestyle='--' , linewidth=1.);
        plt.xticks((0., Sw, Ss, 1.), ['0', '$S_w$', '$S^\star$', '1'])
        plt.yticks((0., 1.), ['0', '1'])
        plt.xlim([0,1+5e-2]); plt.ylim([0,1+5e-2]);
        ax.spines['right'].set_visible(False); ax.spines['top'].set_visible(False) # remove some axes
        ax.plot((1), (0), ls="", marker=">", ms=5, color="k", transform=ax.get_yaxis_transform(), clip_on=False) # make arrows
        ax.plot((0), (1), ls="", marker="^", ms=5, color="k", transform=ax.get_xaxis_transform(), clip_on=False)
        ax.plot([0., Sw, Sw, Ss, Ss, 1.], [0., 0., 0., 1., 1., 1.], color=blue, linewidth=3.0, zorder=10)
        ax.set_aspect('auto', anchor='C')

        # Graph of KR
        ax  = fig.add_subplot(122);
        plt.text(0.38, 0.5, '$K_R(S)$', color=green)
        ax.plot([0., 1.], [1., 1.],  color='k', linestyle='--' , linewidth=1.);
        ax.plot([1., 1.], [0., 1.],  color='k', linestyle='--' , linewidth=1.);
        plt.xticks((0., Sh, 1.), ['0', '$S_h$', '1'])
        plt.yticks((0., 1.), ['0', '1'])
        plt.xlim([0,1+5e-2]); plt.ylim([0,1+5e-2]);
        ax.spines['right'].set_visible(False); ax.spines['top'].set_visible(False) # remove some axes
        ax.plot((1), (0), ls="", marker=">", ms=5, color="k", transform=ax.get_yaxis_transform(), clip_on=False) # make arrows
        ax.plot((0), (1), ls="", marker="^", ms=5, color="k", transform=ax.get_xaxis_transform(), clip_on=False)
        ax.plot([0., Sh, Sh, 1.], [0., 0., 0., 1.], color=green, linewidth=3.0, zorder=10);
        ax.margins(x=5, y=-0.25)

        #plt.savefig("KSKR.png", dpi=120)

wrappers.py

[17]:
!pygmentize wrappers.py
import os
import nutopy as nt
import numpy as np

def singular_control(compile=False, display=False):

    if display:
        out=''
    else:
        out='> /dev/null 2>&1'

    if not compile:
        try:
            from singularcontrol import singularcontrol as sc
        except ModuleNotFoundError:
            compile = True

    if compile:
        # Compiling singular control
        os.system('python -m numpy.f2py -c singularcontrol.f90 \
        h01fun_d01.f90 h0fun_dc0.f90 h1fun_dc1.f90 parameters.f90 \
        -m singularcontrol ' + out)
        from singularcontrol import singularcontrol as sc

    using_split = lambda t, x, p, eta : sc(t, x, p, np.array([eta]))

    @nt.tools.vectorize(vvars=(1, 2, 3))
    def using(t, x, p, eta):
        num, den = using_split(t, x, p, eta)
        return -num/den

    return using

def hamiltonian_umin(compile=False, display=False):

    if display:
        out=''
    else:
        out='> /dev/null 2>&1'

    if not compile:
        try:
            from hmin          import hmin         as hf_min
            from hmin_dmi      import hmin_dmi     as hf_min_d
            from hmin_dmi_dmi  import hmin_dmi_dmi as hf_min_d_d
        except ModuleNotFoundError:
            compile = True

    if compile:
        # Compiling hmin
        os.system('python -m numpy.f2py -c hmin.f90 h0fun.f90 parameters.f90 -m hmin ' + out)
        os.system('python -m numpy.f2py -c hmin_dmi.f90     -m hmin_dmi     ' + out)
        os.system('python -m numpy.f2py -c hmin_dmi_dmi.f90 -m hmin_dmi_dmi ' + out)
        from hmin          import hmin         as hf_min
        from hmin_dmi      import hmin_dmi     as hf_min_d
        from hmin_dmi_dmi  import hmin_dmi_dmi as hf_min_d_d

    # Hamiltonian and flow
    hfun_min   = lambda t, x, p, e                : hf_min(t, x, p, np.array([e]))
    dhfun_min  = lambda t, x, dx, p, dp, e, de    : hf_min_d(t, x, dx, p, dp, np.array([e]), np.array([de]))
    d2hfun_min = lambda t, x, dx, d2x, p, dp, d2p, e, de, d2e : hf_min_d_d(t, x, d2x, dx, p, d2p, dp, \
                                                                           np.array([e]), np.array([d2e]), np.array([de]))

    hfun_min   = nt.tools.tensorize(dhfun_min, d2hfun_min, tvars=(2, 3, 4), full=True)(hfun_min)
    hmin       = nt.ocp.Hamiltonian(hfun_min)

    return hmin

def hamiltonian_umax(compile=False, display=False):

    if display:
        out=''
    else:
        out='> /dev/null 2>&1'

    if not compile:
        try:
            from hmax         import hmax         as hf_max
            from hmax_dma     import hmax_dma     as hf_max_d
            from hmax_dma_dma import hmax_dma_dma as hf_max_d_d
        except ModuleNotFoundError:
            compile = True

    if compile:
        # Compiling hmax
        os.system('python -m numpy.f2py -c hmax.f90 h0fun.f90 h1fun.f90 parameters.f90 -m hmax ' + out)
        os.system('python -m numpy.f2py -c hmax_dma.f90 parameters.f90 -m hmax_dma     ' + out)
        os.system('python -m numpy.f2py -c hmax_dma_dma.f90 hmax_dma.f90 parameters.f90 -m hmax_dma_dma ' + out)
        from hmax         import hmax         as hf_max
        from hmax_dma     import hmax_dma     as hf_max_d
        from hmax_dma_dma import hmax_dma_dma as hf_max_d_d

    # Hamiltonian and flow
    hfun_max   = lambda t, x, p, e                : hf_max(t, x, p, np.array([e]))
    dhfun_max  = lambda t, x, dx, p, dp, e, de    : hf_max_d(t, x, dx, p, dp, np.array([e]), np.array([de]))
    d2hfun_max = lambda t, x, dx, d2x, p, dp, d2p, e, de, d2e : hf_max_d_d(t, x, d2x, dx, p, d2p, dp, \
                                                                           np.array([e]), np.array([d2e]), np.array([de]))

    hfun_max   = nt.tools.tensorize(dhfun_max, d2hfun_max, tvars=(2, 3, 4), full=True)(hfun_max)
    hmax       = nt.ocp.Hamiltonian(hfun_max)

    return hmax

def hamiltonian_using(compile=False, display=False):

    if display:
        out=''
    else:
        out='> /dev/null 2>&1'

    if not compile:
        try:
            from hsing         import hsing         as hf_sing
            from hsing_dus     import hsing_dus     as hf_sing_d
            from hsing_dus_dus import hsing_dus_dus as hf_sing_d_d
        except ModuleNotFoundError:
            compile = True

    if compile:
        # Compiling hsing
        os.system('python -m numpy.f2py -c hsing.f90 h0fun.f90 h1fun.f90 parameters.f90 singularcontrol.f90 \
        h01fun_d01.f90 h0fun_dc0.f90 h1fun_dc1.f90 -m hsing ' + out)
        os.system('python -m numpy.f2py -c hsing_dus.f90 parameters.f90 h1fun_dc1.f90 -m hsing_dus     ' + out)
        os.system('python -m numpy.f2py -c hsing_dus_dus.f90 hsing_dus.f90 parameters.f90 h1fun_dc1.f90 -m hsing_dus_dus ' + out)
        from hsing         import hsing         as hf_sing
        from hsing_dus     import hsing_dus     as hf_sing_d
        from hsing_dus_dus import hsing_dus_dus as hf_sing_d_d

    # Hamiltonian and flow
    hfun_sing   = lambda t, x, p, e                : hf_sing(t, x, p, np.array([e]))
    dhfun_sing  = lambda t, x, dx, p, dp, e, de    : hf_sing_d(t, x, dx, p, dp, np.array([e]), np.array([de]))
    d2hfun_sing = lambda t, x, dx, d2x, p, dp, d2p, e, de, d2e : hf_sing_d_d(t, x, d2x, dx, p, d2p, dp, \
                                                                             np.array([e]), np.array([d2e]), np.array([de]))

    hfun_sing   = nt.tools.tensorize(dhfun_sing, d2hfun_sing, tvars=(2, 3, 4), full=True)(hfun_sing)
    hsing       = nt.ocp.Hamiltonian(hfun_sing)

    return hsing

def switching_function(compile=False, display=False):

    if display:
        out=''
    else:
        out='> /dev/null 2>&1'

    if not compile:
        try:
            from h1fun         import h1fun         as h1f
            from h1fun_dc1     import h1fun_dc1     as h1f_d
        except ModuleNotFoundError:
            compile = True

    if compile:
        # Compiling H1
        os.system('python -m numpy.f2py -c h1fun.f90 parameters.f90 -m h1fun     ' + out)
        os.system('python -m numpy.f2py -c h1fun_dc1.f90 parameters.f90 -m h1fun_dc1 ' + out)
        from h1fun         import h1fun         as h1f
        from h1fun_dc1     import h1fun_dc1     as h1f_d

    h1fun   = lambda p     : h1f(p, np.array([0.0]))
    dh1fun  = lambda p, dp : h1f_d(p, dp, np.array([0.0]))
    h1fun   = nt.tools.tensorize(dh1fun, tvars=(1, ), full=True)(h1fun)
    h1fun   = nt.tools.vectorize(vvars=(1, ))(h1fun)
    #h1      = nt.ocp.Hamiltonian(h1fun)
    h1 = h1fun

    return h1

def lift_H01(compile=False, display=False):

    if display:
        out=''
    else:
        out='> /dev/null 2>&1'

    if not compile:
        try:
            from h01fun         import h01fun         as h01f
            from h01fun_d01     import h01fun_d01     as h01f_d
        except ModuleNotFoundError:
            compile = True

    if compile:
        # Compiling H1
        os.system('python -m numpy.f2py -c h01fun.f90 h0fun_dc0.f90 h1fun_dc1.f90 parameters.f90 -m h01fun     ' + out)
        os.system('python -m numpy.f2py -c h01fun_d01.f90 h1fun_dc1.f90 parameters.f90 -m h01fun_d01 ' + out)
        from h01fun         import h01fun         as h01f
        from h01fun_d01     import h01fun_d01     as h01f_d

    h01fun   = lambda t, x, p, e             : h01f(t, x, p, np.array([e]))
    dh01fun  = lambda t, dt, x, dx, p, dp, e, de : h01f_d(t, dt, x, dx, p, dp, np.array([e]), np.array([de]))
    h01fun   = nt.tools.tensorize(dh01fun, tvars=(1, 2, 3, 4), full=True)(h01fun)
    #h01      = nt.ocp.Hamiltonian(h01fun)
    h01 = h01fun

    return h01

thumbnail