Notebook source code:
examples/irrigation_hybrid/hybrid.ipynb
Run the notebook yourself on binder
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]):
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
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\).
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
Assumption 3. \(k_1\), \(k_2\) are positive parameters with \(k_2 \ge 1\).
Assumption 4. The function \(f\) is the logistic law
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:
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.
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.
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
where the strata are defined by
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)
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
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
Goal 1. We assume that the optimal solution is of the form
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
\(F_1\) being a constant vector field. To solve the optimal control problem by direct shooting, we define the following three different vector fields:
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\):
\(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:
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:
where the cost is given by
where the nonlinear equality constraints are given by
and where the linear inequality constraints are given by
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)
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)
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.
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()
The regularized problem¶
Let \(\eta \ne 0\), \(\eta\) close to \(0\). The regularized optimal control problem, denoted \((\mathrm{OCP}_\eta)\), is the following
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
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:
we have the pseudo-Hamiltonian differential equation
we have the maximization condition
and we have the transversality condition
where
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
However, if \(H_1(t, x(t), p(t)) = 0\) along a time interval on non empty interior, then we have also
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
Derivating twice, we have
Finally, if \(H_1(t, x(t), p(t)) = 0\) along an arc of non empty interior, then the control must satisfy
whenever \(H_{101}(t, z(t)) \ne 0\). Such an arc is called a singular arc. We define the singular control as
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)
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]:
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.
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