Notebook source code:
examples/control-loss/double-integrator-with-control-loss-and-no-control.ipynb
Run the notebook yourself on binder
The double integrator example as a hybrid system¶
Definition of the optimal control problem¶
Here we study a double integrator system as a time optimal control problem. A specifity of our study is that we will consider it as a hybrid system, namely a dynamic system which dynamics changes depending on some conditions. In our case, the system goes through three stages:
In the first region (\(\mathbb{R}⁻\setminus[-0.65, 0]\)), we lose the ability to change the control of our system: the control becomes a constant, and has to be in [-1, 1];
In the second region (\([-0.65, 0]\)), the control is set to zero and cannot be changed;
In the third region, we can change the control at will, and has to be in [-1, 1].
The corresponding problem is defined below:
In the following, we will consider \(\lambda\) as a constant parameter of the system. In order to do that, we will take it as another state of the system, with no dynamics.
General idea¶
Here we offer a method to solve this kind of problem by using first direct solving and then indirect solving, which is the standard. The main difference is that we also have to do additional analysis regarding the different phases before doing indirect solving.
[1]:
# Imports for direct methods
using JuMP # NLP modeling
using Ipopt # NLP solving
# Imports to plot solutions
using Plots
using Plots.PlotMeasures
[2]:
mutable struct DirectSolution
t; x1; x2; x3; p1; p2 ; u # x3 denotes the extended state i.e. λ
end;
[3]:
# Parameters
t0 = 0. # initial time
x1f = 0. # Final position
x2f = 0.
M = 1.0 # control bound
no_control_bound = -0.65; # lower bound for no control region
[4]:
# Heaviside regularization
function G(x, threshold=0.)
n = 10. # regularization stiffness
if - 1/(2*n) < x - threshold < 1/(2*n)
return 1. - 3*abs2(n*(x-threshold) + 0.5) + 2*abs2(n*(x-threshold) + 0.5)*(n*(x-threshold) + 0.5)
elseif x - threshold ≤ - 1/(2*n)
return 1.
else x - threshold ≥ 1/(2*n)
return 0.
end
end;
[5]:
function DI(x0; solution=[], nsteps=200, display=true)
# Create JuMP model, using Ipopt as the solver
if display
pl = 5
else
pl = 1
end
sys = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => pl))
set_optimizer_attribute(sys,"tol",1e-8)
set_optimizer_attribute(sys,"constr_viol_tol",1e-8)
set_optimizer_attribute(sys,"max_iter",1000)
# register G function with 1 argument
register(sys, :G, 2, G; autodiff = true)
N = nsteps # Grid size
@variables(sys, begin
x1[1:N+1] # x1
x2[1:N+1] # x2
x3[1:N+1] ##### x3=λ #####
-M ≤ u[i=1:N+1] ≤ M # Control
0. ≤ Δt ≤ 1. # contraints may be "superflues"
end)
# Objective
@objective(sys, Min, Δt)
# Boundary constraints
@constraints(sys, begin
con_x10, x1[1] - x0[1] == 0
con_x20, x2[1] - x0[2] == 0
con_x30, -M ≤ x3[1] ≤ M # parameter lambda constraint
con_x1f, x1[N+1] - x1f == 0
con_x2f, x2[N+1] - x2f == 0
con_x3f, -M ≤ x3[N+1] ≤ M # parameter lambda constraint
end)
# Dynamics
@NLexpression(sys, dx1[j = 1:N+1], x2[j])
@NLexpression(sys, dx2[j = 1:N+1], x3[j]*G(x1[j], no_control_bound) + (1. - G(x1[j], 0.))*u[j]) # regularized dynamics
@NLexpression(sys, dx3[j = 1:N+1], 0.) # augmented dynamics
# Dynamics with Crank-Nicolson scheme
@NLconstraints(sys, begin
con_dx1[j=1:N], x1[j+1] == x1[j] + 0.5 * Δt * (dx1[j+1] + dx1[j])
con_dx2[j=1:N], x2[j+1] == x2[j] + 0.5 * Δt * (dx2[j+1] + dx2[j])
con_dx3[j=1:N], x3[j+1] == x3[j]
end);
# Solve for the control and state
if display
println("Solving...")
end
status = optimize!(sys)
if display
println()
end
# Display results
if display
if termination_status(sys) == MOI.OPTIMAL
println(" Solution is optimal")
elseif termination_status(sys) == MOI.LOCALLY_SOLVED
println(" (Local) solution found")
elseif termination_status(sys) == MOI.TIME_LIMIT && has_values(sys)
println(" Solution is suboptimal due to a time limit, but a primal solution is available")
else
error(" The model was not solved correctly.")
end
println(" objective value = ", objective_value(sys))
println()
end
# Retrieves values (including duals)
x1 = value.(x1)[:]
x2 = value.(x2)[:]
x3 = value.(x3)[:]
u = value.(u)[:]
t = (0:N) * value.(Δt)
px10 = dual(con_x10)
px20 = dual(con_x20)
px1f = dual(con_x1f)
px2f = dual(con_x2f)
if(px10*dual(con_dx1[1])<0); px10 = -px10; end
if(px20*dual(con_dx2[1])<0); px20 = -px20; end
if(px1f*dual(con_dx1[N])<0); px1f = -px1f; end
if(px2f*dual(con_dx2[N])<0); px2f = -px2f; end
# H = px1 * x2 + px2 * u > 0, according to Pontryagin Maximum Principle convention
if (px10 * x2[1] + px20 * u[1] < 0)
sign = -1.0 # switch sign of dual variables
else
sign = 1.0
end
px1 = [ dual(con_dx1[i]) for i in 1:N ]
px2 = [ dual(con_dx2[i]) for i in 1:N ]
px1 = sign * [px10; (px1[1:N-1]+px1[2:N])/2; px1f]; # We add the multiplier from the limit conditions
px2 = sign * [px20; (px2[1:N-1]+px2[2:N])/2; px2f]; # We add the multiplier from the limit conditions
p1 = sign * [px10; px1[1:N-1]; px1f]; # We add the multiplier from the limit conditions
p2 = sign * [px20; px2[1:N-1]; px2f]; # We add the multiplier from the limit conditions
return DirectSolution(t, x1, x2, x3, p1, p2, u)
end;
[6]:
# Solving
x0_ref = [-1.; 0.; ]
direct_sol= DI(x0_ref);
Solving...
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit https://github.com/coin-or/Ipopt
******************************************************************************
This is Ipopt version 3.14.4, running with linear solver MUMPS 5.4.1.
Number of nonzeros in equality constraint Jacobian...: 3204
Number of nonzeros in inequality constraint Jacobian.: 4
Number of nonzeros in Lagrangian Hessian.............: 0
Total number of variables............................: 805
variables with only lower bounds: 0
variables with lower and upper bounds: 202
variables with only upper bounds: 0
Total number of equality constraints.................: 604
Total number of inequality constraints...............: 4
inequality constraints with only lower bounds: 2
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 2
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 9.9999900e-03 1.00e+00 1.00e+00 0.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 6.8099940e-03 6.75e-01 2.55e+04 -6.0 3.04e+00 - 2.45e-01 3.25e-01h 1
2 2.0725749e-02 3.66e-01 6.03e+04 0.8 3.97e+00 - 1.00e+00 4.57e-01f 1
3 3.6418072e-02 1.90e-01 2.30e+05 0.7 3.66e+00 - 1.03e-01 4.81e-01f 1
4 3.3868178e-02 2.79e-02 9.59e+02 0.7 1.15e+00 - 1.00e+00 1.00e+00f 1
5 4.1000889e-02 2.04e-02 7.30e+03 1.3 3.06e+00 - 4.02e-01 2.71e-01h 1
6 5.4020703e-02 2.12e-02 5.10e+00 0.6 2.66e-01 - 9.99e-01 1.00e+00f 1
7 8.0224667e-02 8.48e-03 4.60e+01 0.5 8.67e-01 - 1.00e+00 1.00e+00f 1
8 1.4124039e-01 3.41e-03 4.20e+00 0.5 1.27e-01 - 9.25e-01 1.00e+00f 1
9 2.4323902e-01 4.38e-03 2.56e+00 -0.2 1.02e-01 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 2.0595934e-01 8.34e-04 1.77e-01 -0.9 3.73e-02 - 1.00e+00 1.00e+00h 1
11 1.7049810e-01 6.60e-03 3.71e+00 -6.9 9.11e-02 - 9.23e-01 1.00e+00h 1
12 1.0638898e-01 4.09e-03 8.60e-01 -3.4 6.41e-02 - 1.00e+00 1.00e+00h 1
13 1.0195017e-01 1.96e-04 3.44e-01 -4.7 7.04e-02 - 1.00e+00 1.00e+00h 1
14 9.6217689e-02 8.72e-05 7.76e-02 -6.3 2.08e-02 - 1.00e+00 1.00e+00h 1
15 6.2270123e-02 3.17e-03 2.62e-01 -7.3 1.27e-01 - 1.00e+00 1.00e+00h 1
16 4.1001332e-02 7.14e-03 2.87e-02 -8.2 3.44e-01 - 1.00e+00 1.00e+00h 1
17 4.1948638e-02 6.98e-04 2.12e-01 -6.0 3.85e-01 - 1.00e+00 9.28e-01h 1
18 4.1290979e-02 5.90e-04 2.02e-01 -4.8 1.90e-01 - 1.00e+00 1.45e-01h 1
19 3.9319057e-02 3.77e-04 1.72e-01 -4.5 4.98e-01 - 1.00e+00 2.18e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20 3.7159893e-02 5.93e-04 1.96e-01 -4.8 9.39e-01 - 1.00e+00 2.02e-01h 1
21 3.4899661e-02 1.81e-03 3.20e-01 -5.1 1.99e+00 - 1.00e+00 1.65e-01h 1
22 3.2125635e-02 6.31e-03 7.48e-01 -4.8 5.55e+00 - 1.00e+00 1.23e-01h 1
23 3.1520033e-02 4.29e-04 6.43e-01 -4.3 3.96e-01 - 1.00e+00 1.00e+00h 1
24 3.0367934e-02 1.16e-03 5.67e-01 -4.7 5.04e-01 - 1.00e+00 8.41e-01h 1
25 2.9959235e-02 1.29e-03 4.14e-01 -4.4 6.16e-01 - 1.00e+00 4.85e-01h 1
26 2.9847837e-02 2.54e-06 1.81e-02 -5.7 2.55e-03 - 1.00e+00 1.00e+00h 1
27 2.9509466e-02 8.66e-06 4.17e-03 -6.7 1.06e-01 - 1.00e+00 2.95e-01h 1
28 2.9253244e-02 1.27e-05 2.56e-03 -7.1 8.76e-01 - 1.00e+00 7.73e-02h 1
29 2.9055736e-02 6.11e-05 1.71e-03 -7.3 7.98e-01 - 1.00e+00 1.38e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
30 2.8846682e-02 1.41e-04 1.87e-03 -7.1 1.99e+00 - 1.00e+00 1.05e-01h 1
31 2.8734906e-02 1.60e-04 1.53e-03 -5.9 2.11e+00 - 1.00e+00 7.31e-02h 1
32 2.8558178e-02 2.32e-04 1.38e-03 -5.4 4.16e+00 - 1.00e+00 7.23e-02h 1
33 2.8217358e-02 8.99e-04 4.85e-03 -5.2 2.70e+00 - 1.00e+00 2.69e-01h 1
34 2.7984517e-02 1.53e-03 6.46e-03 -4.7 3.94e+00 - 1.00e+00 2.27e-01h 1
35 2.7864914e-02 3.47e-03 3.38e-02 -4.4 2.90e+00 - 1.00e+00 3.79e-01h 1
36 2.7532957e-02 3.57e-03 7.41e-02 -4.4 2.23e+00 - 9.99e-01 4.31e-01h 1
37 2.7328035e-02 2.00e-05 7.12e-02 -5.5 8.11e-03 - 1.00e+00 1.00e+00h 1
38 2.6846330e-02 2.30e-05 4.33e-02 -6.8 7.40e-02 - 1.00e+00 6.99e-01h 1
39 2.6709988e-02 2.91e-05 4.88e-02 -7.3 9.78e-01 - 1.00e+00 2.00e-02h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
40 2.6310453e-02 3.89e-04 9.53e-02 -7.4 2.58e+00 - 1.00e+00 7.61e-02h 1
41 2.6587368e-02 2.02e-03 8.41e-02 -4.6 9.11e-01 - 1.00e+00 6.91e-01h 1
42 2.6413207e-02 9.83e-04 1.56e-01 -4.7 9.04e-01 - 9.17e-01 9.20e-01h 1
43 2.6115710e-02 5.56e-03 2.82e-01 -4.5 1.26e+00 - 1.00e+00 5.71e-01f 1
44 2.5947214e-02 3.57e-03 2.54e-01 -4.2 1.27e+00 - 1.00e+00 3.55e-01h 1
45 2.5711710e-02 2.76e-05 1.68e-01 -5.4 7.80e-03 - 1.00e+00 1.00e+00h 1
46 2.5284666e-02 5.71e-05 5.98e-02 -6.8 4.95e-02 - 1.00e+00 1.00e+00h 1
47 2.4605960e-02 1.07e-04 6.07e-02 -7.4 1.22e+00 - 1.00e+00 8.03e-02h 1
48 2.4319683e-02 1.60e-04 1.91e-02 -6.4 1.06e+00 - 1.00e+00 7.23e-02h 1
49 2.4061624e-02 3.37e-04 4.51e-02 -6.2 1.30e+00 - 1.00e+00 1.13e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
50 2.3886736e-02 1.16e-03 1.67e-01 -5.1 9.93e-01 - 1.00e+00 3.93e-01h 1
51 2.3844639e-02 7.08e-04 1.92e-01 -4.9 8.38e-01 - 1.00e+00 5.13e-01h 1
52 2.3820378e-02 1.11e-06 6.31e-03 -6.0 1.57e-03 - 1.00e+00 1.00e+00h 1
53 2.3662491e-02 1.36e-05 1.14e-02 -6.5 5.20e-01 - 1.00e+00 9.54e-02h 1
54 2.3471776e-02 3.77e-05 1.02e-02 -7.6 4.83e-01 - 1.00e+00 3.38e-01h 1
55 2.3347765e-02 1.46e-04 4.91e-03 -6.3 1.21e+00 - 1.00e+00 1.90e-01h 1
56 2.3294937e-02 1.64e-04 8.34e-03 -6.0 2.18e+00 - 1.00e+00 5.64e-02h 1
57 2.3192120e-02 2.48e-04 1.15e-02 -5.7 3.47e+00 - 1.00e+00 1.01e-01h 1
58 2.3106748e-02 3.52e-04 2.02e-03 -5.5 2.60e+00 - 1.00e+00 1.72e-01h 1
59 2.3039015e-02 6.14e-04 1.33e-02 -5.0 4.65e+00 - 1.00e+00 1.33e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
60 2.3465032e-02 3.36e-03 8.16e-02 -4.6 2.03e+00 - 1.00e+00 6.70e-01h 1
61 2.3253455e-02 1.71e-03 2.69e-02 -4.7 9.62e-01 - 9.44e-01 7.82e-01h 1
62 2.3087870e-02 1.44e-03 3.71e-02 -4.7 1.07e+00 - 1.00e+00 4.91e-01f 1
63 2.2885272e-02 1.78e-03 4.13e-02 -4.7 6.79e-01 - 1.00e+00 7.06e-01f 1
64 2.2841142e-02 4.17e-06 2.34e-02 -5.9 8.40e-03 - 1.00e+00 1.00e+00h 1
65 2.2572466e-02 4.05e-05 6.17e-02 -6.6 3.95e-01 - 1.00e+00 1.04e-01h 1
66 2.2392367e-02 3.34e-05 3.36e-02 -6.5 3.25e-01 - 1.00e+00 1.81e-01h 1
67 2.2131111e-02 4.37e-05 6.26e-02 -6.5 1.64e+00 - 1.00e+00 9.29e-02h 1
68 2.2004578e-02 4.65e-05 7.64e-02 -6.0 1.77e+00 - 1.00e+00 2.49e-01h 1
69 2.1964143e-02 4.10e-05 6.34e-02 -6.1 1.18e+00 - 1.00e+00 1.62e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
70 2.1895274e-02 8.05e-05 2.41e-02 -6.1 1.33e+00 - 1.00e+00 4.68e-01h 1
71 2.1867797e-02 1.01e-04 1.62e-02 -6.1 3.09e+00 - 1.00e+00 9.54e-02h 1
72 2.1713170e-02 4.83e-04 3.49e-02 -6.1 1.01e+00 - 1.00e+00 9.19e-01h 1
73 2.1675720e-02 4.48e-04 3.39e-02 -11.0 2.68e+00 - 3.58e-01 1.69e-01h 1
74 2.1528096e-02 3.47e-04 3.53e-02 -6.4 6.72e-01 - 1.00e+00 9.98e-01h 1
75 2.1531496e-02 2.78e-08 1.82e-04 -7.6 1.13e-03 - 1.00e+00 1.00e+00h 1
76 2.1529603e-02 4.87e-10 2.39e-04 -9.6 2.40e-04 - 1.00e+00 1.00e+00h 1
77 2.1512198e-02 1.46e-06 1.65e-03 -8.7 4.10e-01 - 1.00e+00 8.48e-03f 1
78 2.1498874e-02 4.63e-06 2.52e-03 -8.1 4.33e-01 - 1.00e+00 2.91e-02h 1
79 2.1476301e-02 2.33e-05 7.27e-03 -6.9 2.05e+00 - 1.00e+00 2.20e-02h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
80 2.1378827e-02 5.95e-04 4.00e-02 -5.1 2.76e+01 - 1.00e+00 1.43e-02h 1
81 2.1269391e-02 1.27e-03 3.10e-02 -7.0 3.87e+00 - 2.36e-01 2.03e-01h 1
82 2.1545582e-02 2.09e-03 1.51e-01 -5.0 2.76e+00 - 9.63e-01 6.57e-01h 1
83 2.1546382e-02 3.80e-07 4.75e-03 -5.1 4.75e-03 - 1.00e+00 1.00e+00h 1
84 2.1459461e-02 3.33e-04 2.30e-02 -5.1 2.04e-01 - 1.00e+00 1.00e+00h 1
85 2.1767337e-02 1.08e-04 3.60e-03 -4.8 3.69e-01 - 1.00e+00 7.44e-01h 1
86 2.1673182e-02 3.66e-04 7.80e-03 -4.9 6.62e-01 - 8.35e-01 8.94e-01h 1
87 2.1585159e-02 7.01e-04 8.83e-03 -4.9 1.05e+00 - 1.00e+00 8.30e-01f 1
88 2.1487690e-02 1.09e-03 1.15e-02 -4.9 8.09e-01 - 1.00e+00 8.47e-01f 1
89 2.1490895e-02 1.88e-06 2.67e-01 -4.9 2.67e-01 - 1.00e+00 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
90 2.1371133e-02 2.19e-05 8.89e-03 -4.9 1.22e+00 - 9.90e-01 3.66e-01h 1
91 2.0774678e-02 1.55e-04 5.47e-02 -5.6 5.63e-01 - 1.00e+00 1.00e+00h 1
92 2.0669360e-02 3.74e-05 1.55e-02 -5.9 2.40e-01 - 1.00e+00 9.90e-01h 1
93 2.0580958e-02 1.10e-05 1.36e-02 -6.7 1.91e-01 - 1.00e+00 1.00e+00h 1
94 2.0551307e-02 5.32e-05 3.93e-02 -7.1 1.33e+00 - 1.00e+00 2.32e-01h 1
95 2.0527457e-02 1.19e-04 6.39e-02 -6.2 1.20e+00 - 1.00e+00 3.19e-01h 1
96 2.0509784e-02 1.03e-04 6.42e-02 -11.0 3.17e-01 - 6.25e-01 1.93e-01h 1
97 2.0507903e-02 9.51e-09 1.72e-04 -8.1 1.79e-04 - 1.00e+00 1.00e+00h 1
98 2.0492402e-02 4.42e-06 4.30e-03 -7.0 1.67e-01 - 1.00e+00 3.25e-02h 1
99 2.0297403e-02 2.85e-03 1.26e-02 -6.1 5.32e-01 - 1.00e+00 8.37e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
100 2.0455861e-02 2.14e-04 2.94e-02 -5.7 1.90e-01 - 1.00e+00 1.00e+00h 1
101 2.0376003e-02 1.27e-04 9.22e-03 -6.1 3.29e-01 - 1.00e+00 1.00e+00h 1
102 2.0275631e-02 2.71e-04 1.78e-02 -7.1 3.77e-01 - 1.00e+00 9.99e-01h 1
103 2.0260456e-02 1.40e-04 1.35e-03 -6.3 1.54e-01 - 1.00e+00 1.00e+00h 1
104 2.0260023e-02 3.58e-09 1.88e-04 -8.4 3.21e-04 - 1.00e+00 1.00e+00h 1
105 2.0242758e-02 6.54e-07 2.38e-03 -9.0 4.32e-01 - 1.00e+00 2.49e-02h 1
106 2.0232234e-02 5.56e-06 7.25e-03 -9.0 2.01e+00 - 1.00e+00 3.77e-02h 1
107 2.0202825e-02 5.28e-05 2.23e-02 -8.5 1.96e+00 - 1.00e+00 1.18e-01h 1
108 2.0193367e-02 5.68e-05 2.67e-02 -7.2 1.76e+00 - 1.00e+00 3.71e-02h 1
109 2.0191181e-02 9.88e-09 1.85e-04 -8.5 1.12e-04 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
110 2.0189715e-02 7.54e-07 1.27e-03 -6.8 9.64e-02 - 1.00e+00 7.03e-02h 1
111 2.0088908e-02 2.51e-03 3.13e-02 -6.6 9.28e-01 - 1.00e+00 7.28e-01h 1
112 2.0308218e-02 5.82e-03 4.57e-02 -5.6 3.62e+00 - 1.00e+00 4.56e-01h 1
113 2.0322750e-02 1.48e-03 3.03e-02 -5.6 7.09e-01 - 1.00e+00 1.00e+00h 1
114 2.0219050e-02 2.37e-03 7.62e-02 -5.6 1.25e+00 - 6.17e-01 4.19e-01h 1
115 2.0188692e-02 1.87e-05 1.57e-01 -5.6 1.57e-01 - 1.00e+00 1.00e+00f 1
116 2.0189051e-02 4.94e-10 1.89e-04 -5.6 2.11e-04 - 1.00e+00 1.00e+00h 1
117 2.0187531e-02 5.69e-09 2.35e-04 -11.0 2.36e-04 - 1.00e+00 1.00e+00h 1
118 2.0185992e-02 5.85e-09 2.41e-04 -11.0 2.41e-04 - 1.00e+00 1.00e+00h 1
119 2.0132979e-02 4.59e-07 3.38e-03 -10.3 1.36e+00 - 1.00e+00 6.35e-03f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
120 2.0085731e-02 7.15e-07 3.04e-03 -10.7 6.81e-01 - 1.00e+00 1.63e-02h 1
121 2.0049577e-02 4.44e-06 5.12e-03 -11.0 1.82e+00 - 1.00e+00 1.66e-02h 1
122 2.0010279e-02 3.46e-05 1.34e-02 -9.7 1.37e+01 - 1.00e+00 6.11e-03h 1
123 1.9970532e-02 9.16e-05 2.53e-02 -7.7 9.87e+00 - 1.00e+00 1.47e-02h 1
124 1.9936853e-02 1.22e-04 3.55e-02 -7.1 6.92e+00 - 1.00e+00 2.15e-02h 1
125 1.9935313e-02 2.31e-08 2.24e-04 -7.9 1.31e-04 - 1.00e+00 1.00e+00h 1
126 1.9924046e-02 1.49e-07 1.15e-03 -9.0 1.46e-01 - 1.00e+00 5.61e-02h 1
127 1.9913805e-02 2.73e-07 2.96e-03 -8.9 2.00e+00 - 1.00e+00 2.92e-02h 1
128 1.9877966e-02 7.79e-06 7.13e-03 -7.3 8.43e+01 - 1.00e+00 4.42e-03h 1
129 1.9847507e-02 1.47e-05 1.23e-02 -6.5 1.49e+01 - 1.00e+00 2.19e-02h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
130 1.9820287e-02 2.58e-05 2.13e-02 -6.9 6.39e+00 - 1.00e+00 1.25e-01h 1
131 1.9819131e-02 2.21e-09 1.48e-04 -7.6 1.04e-04 - 1.00e+00 1.00e+00h 1
132 1.9809282e-02 6.81e-05 1.27e-02 -6.3 2.66e-01 - 9.98e-01 2.05e-01h 1
133 1.9783323e-02 5.58e-04 9.34e-03 -6.1 2.57e-01 - 1.00e+00 5.96e-01h 1
134 1.9762729e-02 4.99e-05 1.25e-02 -6.8 1.78e-01 - 1.00e+00 1.00e+00h 1
135 1.9735015e-02 8.17e-05 1.16e-02 -7.7 2.13e-01 - 1.00e+00 1.00e+00h 1
136 1.9707892e-02 1.52e-04 1.38e-02 -7.8 2.58e-01 - 1.00e+00 1.00e+00h 1
137 1.9705549e-02 3.01e-08 1.37e-04 -8.0 2.49e-04 - 1.00e+00 1.00e+00h 1
138 1.9705319e-02 2.97e-08 8.17e-05 -9.3 8.58e-02 - 1.00e+00 1.43e-02h 1
139 1.9694991e-02 3.05e-07 1.45e-03 -9.9 8.53e-02 - 1.00e+00 8.95e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
140 1.9687512e-02 3.99e-06 2.29e-03 -8.4 8.60e-02 - 1.00e+00 7.17e-01h 1
141 1.9687257e-02 4.03e-11 7.80e-05 -9.6 7.60e-05 - 1.00e+00 1.00e+00h 1
142 1.9679185e-02 9.40e-06 2.70e-03 -10.2 4.77e-01 - 1.00e+00 1.38e-01h 1
143 1.9663956e-02 3.95e-05 7.21e-03 -8.2 4.97e-01 - 1.00e+00 2.38e-01h 1
144 1.9643317e-02 1.94e-04 1.99e-02 -7.5 2.53e-01 - 1.00e+00 1.00e+00h 1
145 1.9641542e-02 2.66e-08 1.23e-04 -8.0 3.53e-04 - 1.00e+00 1.00e+00h 1
146 1.9639812e-02 1.32e-07 1.04e-04 -9.2 1.21e-01 - 1.00e+00 5.68e-02h 1
147 1.9637820e-02 4.05e-07 5.38e-04 -10.3 7.90e-02 - 1.00e+00 1.33e-01h 1
148 1.9628368e-02 1.29e-05 3.41e-03 -9.3 8.02e-02 - 1.00e+00 9.83e-01h 1
149 1.9628078e-02 1.14e-10 6.80e-05 -10.8 6.38e-05 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
150 1.9624731e-02 1.97e-06 1.30e-03 -11.0 3.00e-01 - 1.00e+00 1.02e-01h 1
151 1.9600292e-02 1.11e-04 1.07e-02 -11.0 3.07e-01 - 1.00e+00 7.16e-01h 1
152 1.9600251e-02 1.11e-04 1.07e-02 -8.8 3.35e-01 - 1.00e+00 8.45e-04h 1
153 1.9600119e-02 1.66e-08 3.74e-03 -7.6 3.70e-03 - 1.00e+00 1.00e+00h 1
154 1.9569130e-02 3.13e-04 1.65e-02 -7.7 1.02e+00 - 1.00e+00 4.08e-01h 1
155 1.9570485e-02 4.06e-04 2.21e-02 -5.7 1.31e+00 - 1.00e+00 1.95e-01h 1
156 1.9591186e-02 6.11e-04 2.17e-02 -5.8 7.93e-01 - 1.00e+00 8.70e-01h 1
157 1.9594396e-02 3.15e-06 2.38e-01 -5.8 2.38e-01 - 1.00e+00 1.00e+00f 1
158 1.9597537e-02 9.46e-06 2.28e-03 -5.4 2.61e+00 - 9.90e-01 1.59e-01h 1
159 1.9698976e-02 2.10e-05 5.78e-03 -5.4 5.55e-01 - 1.00e+00 6.99e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
160 1.9742880e-02 1.49e-05 1.49e-02 -5.4 1.58e-01 - 1.00e+00 1.00e+00f 1
161 1.9741124e-02 1.77e-06 1.60e-02 -5.4 3.14e-02 - 1.00e+00 1.00e+00h 1
162 1.9738438e-02 1.71e-06 1.60e-02 -5.4 4.33e-02 - 1.00e+00 1.00e+00h 1
163 1.9738466e-02 1.13e-09 2.36e-05 -5.4 2.42e-05 - 1.00e+00 1.00e+00h 1
164 1.9513439e-02 6.81e-04 2.39e-02 -8.1 9.54e-02 - 9.17e-01 7.36e-01f 1
165 1.9751229e-02 4.35e-04 1.33e-02 -5.3 1.96e-01 - 1.00e+00 1.00e+00f 1
166 1.9687258e-02 4.05e-04 1.36e-02 -5.4 6.39e-01 - 1.00e+00 4.74e-01H 1
167 1.9605931e-02 2.52e-04 1.35e-02 -5.4 5.90e-01 - 1.00e+00 1.00e+00f 1
168 1.9592304e-02 2.50e-04 1.41e-02 -5.4 4.56e-01 - 1.00e+00 2.05e-01h 2
169 1.9587950e-02 1.32e-07 1.98e-04 -6.9 3.46e-04 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
170 1.9493741e-02 3.92e-06 6.98e-03 -7.9 1.92e-01 - 1.00e+00 7.07e-02h 1
171 1.9393816e-02 2.31e-05 1.04e-02 -8.4 6.10e-01 - 1.00e+00 9.86e-02h 1
172 1.9372381e-02 4.41e-05 1.59e-02 -8.5 2.38e+00 - 1.00e+00 2.50e-02h 1
173 1.9335741e-02 2.71e-04 3.01e-02 -6.2 1.39e+00 - 1.00e+00 2.25e-01h 1
174 1.9788526e-02 8.94e-04 7.41e-02 -5.0 9.59e-01 - 1.00e+00 7.51e-01h 1
175 1.9796739e-02 7.46e-08 2.72e-03 -5.1 2.71e-03 - 1.00e+00 1.00e+00h 1
176 1.9796736e-02 3.45e-10 8.42e-04 -5.1 8.42e-04 - 1.00e+00 1.00e+00h 1
177 1.9796737e-02 3.46e-10 9.12e-04 -5.1 9.12e-04 - 1.00e+00 1.00e+00h 1
178 1.9796729e-02 3.28e-10 7.93e-04 -5.1 7.93e-04 - 1.00e+00 1.00e+00h 1
179 1.9796717e-02 3.13e-10 7.24e-04 -5.1 7.24e-04 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
180 1.9796700e-02 2.98e-10 6.67e-04 -5.1 6.67e-04 - 1.00e+00 1.00e+00h 1
181 1.9796680e-02 2.85e-10 6.20e-04 -5.1 6.20e-04 - 1.00e+00 1.00e+00h 1
182 1.9796657e-02 2.73e-10 5.81e-04 -5.1 5.81e-04 - 1.00e+00 1.00e+00h 1
183 1.9796632e-02 2.61e-10 5.48e-04 -5.1 5.48e-04 - 1.00e+00 1.00e+00h 1
184 1.9796604e-02 2.51e-10 5.19e-04 -5.1 5.19e-04 - 1.00e+00 1.00e+00h 1
185 1.9796576e-02 2.41e-10 4.93e-04 -5.1 4.93e-04 - 1.00e+00 1.00e+00h 1
186 1.9796546e-02 2.32e-10 4.71e-04 -5.1 4.71e-04 - 1.00e+00 1.00e+00h 1
187 1.9796514e-02 2.24e-10 4.51e-04 -5.1 4.51e-04 - 1.00e+00 1.00e+00h 1
188 1.9796482e-02 2.17e-10 4.33e-04 -5.1 4.33e-04 - 1.00e+00 1.00e+00h 1
189 1.9796448e-02 2.10e-10 4.17e-04 -5.1 4.17e-04 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
190 1.9796414e-02 2.03e-10 4.02e-04 -5.1 4.02e-04 - 1.00e+00 1.00e+00h 1
191 1.9796379e-02 1.97e-10 3.89e-04 -5.1 3.89e-04 - 1.00e+00 1.00e+00h 1
192 1.9762208e-02 5.71e-04 3.33e-02 -5.1 1.25e+00 - 1.00e+00 2.34e-01f 2
193 1.9739419e-02 1.43e-04 5.11e-02 -5.1 2.77e-01 - 1.00e+00 1.00e+00h 1
194 1.9726498e-02 3.03e-04 6.01e-02 -5.1 2.35e-01 - 1.00e+00 1.00e+00h 1
195 1.9681855e-02 1.39e-03 8.29e-02 -5.1 5.87e-01 - 1.00e+00 1.00e+00h 1
196 1.9657173e-02 3.20e-06 3.15e-02 -5.1 2.72e-03 - 1.00e+00 1.00e+00h 1
197 1.9657062e-02 3.80e-10 1.35e-04 -5.1 1.36e-04 - 1.00e+00 1.00e+00h 1
198 1.9627775e-02 5.98e-06 2.46e-03 -5.1 1.04e-01 - 1.00e+00 4.56e-01h 2
199 1.9613721e-02 5.24e-06 2.28e-03 -5.1 2.62e-01 - 1.00e+00 4.67e-01h 2
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
200 1.9603956e-02 2.73e-05 3.48e-03 -5.1 1.48e-01 - 1.00e+00 1.00e+00h 1
201 1.9598604e-02 2.10e-05 5.04e-03 -5.1 1.22e-01 - 1.00e+00 1.00e+00h 1
202 1.9598424e-02 4.39e-11 1.46e-05 -5.1 3.86e-05 - 1.00e+00 1.00e+00h 1
203 1.9597280e-02 4.88e-09 1.89e-04 -11.0 1.90e-04 - 1.00e+00 1.00e+00h 1
204 1.9434969e-02 8.87e-06 8.45e-03 -7.0 1.15e+00 - 1.00e+00 2.21e-02f 1
205 1.9265961e-02 2.68e-05 6.35e-03 -6.7 6.16e-01 - 1.00e+00 7.92e-02h 1
206 1.9201749e-02 6.17e-05 1.26e-02 -7.0 6.42e-01 - 1.00e+00 6.71e-02h 1
207 1.9158793e-02 1.62e-04 1.85e-02 -6.6 1.55e+00 - 1.00e+00 9.72e-02h 1
208 1.9129191e-02 7.76e-04 1.83e-02 -5.8 9.40e-01 - 1.00e+00 6.81e-01h 1
209 1.9238446e-02 1.51e-05 6.36e-03 -5.6 7.04e-01 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
210 1.9238444e-02 4.39e-09 1.55e-04 -6.4 1.55e-04 - 1.00e+00 1.00e+00h 1
211 1.9097276e-02 3.10e-06 9.62e-03 -7.3 1.51e-01 - 1.00e+00 1.57e-01f 1
212 1.9038135e-02 9.39e-05 1.13e-02 -7.0 6.99e-01 - 1.00e+00 3.39e-01h 1
213 1.9028909e-02 5.77e-04 1.84e-02 -6.0 7.07e-01 - 1.00e+00 6.30e-01h 1
214 1.8992378e-02 1.24e-04 1.52e-02 -6.4 1.93e-01 - 1.00e+00 1.00e+00h 1
215 1.8992882e-02 1.66e-09 1.35e-04 -7.6 2.57e-04 - 1.00e+00 1.00e+00h 1
216 1.8972075e-02 3.87e-07 1.56e-03 -7.9 1.01e+00 - 1.00e+00 4.67e-02h 1
217 1.8957577e-02 1.10e-06 1.86e-03 -8.3 2.17e+00 - 1.00e+00 1.46e-01h 1
218 1.8946424e-02 1.56e-06 3.32e-03 -7.8 4.91e+01 - 1.00e+00 6.91e-03h 1
219 1.8939970e-02 2.02e-06 5.00e-03 -7.0 7.68e+00 - 1.00e+00 1.11e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
220 1.8928607e-02 1.71e-06 7.88e-03 -7.8 2.84e-01 - 1.00e+00 7.57e-01h 1
221 1.8928032e-02 9.41e-10 8.78e-05 -8.7 8.75e-05 - 1.00e+00 1.00e+00h 1
222 1.8926329e-02 7.86e-06 3.23e-03 -7.0 2.84e-01 - 1.00e+00 6.37e-02h 1
223 1.8870014e-02 7.01e-03 1.04e-02 -6.1 6.11e-01 - 8.19e-01 1.00e+00f 1
224 1.8673177e-02 4.87e-03 6.02e-03 -6.9 3.84e-01 - 1.00e+00 1.00e+00h 1
225 1.8917200e-02 1.26e-04 4.43e-03 -6.5 3.96e-02 - 1.00e+00 1.00e+00h 1
226 1.8911722e-02 6.14e-05 1.95e-03 -6.5 1.73e-01 - 1.00e+00 9.91e-01h 1
227 1.8904890e-02 1.60e-03 4.32e-02 -5.1 3.09e+01 - 1.00e+00 2.43e-02h 1
228 2.0014284e-02 1.13e-03 4.94e-02 -4.7 1.52e+00 - 1.00e+00 8.04e-01h 1
229 1.9805486e-02 8.82e-04 5.86e-02 -4.9 5.02e-01 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
230 1.9803428e-02 5.42e-06 6.86e-02 -4.9 2.34e-01 - 1.00e+00 1.00e+00h 1
231 1.9801188e-02 2.06e-04 4.55e-02 -4.9 1.45e-01 - 1.00e+00 1.00e+00h 1
232 1.9793936e-02 6.83e-05 2.69e-02 -4.9 2.92e-01 - 1.00e+00 1.00e+00h 1
233 1.9792660e-02 1.48e-08 2.39e-05 -4.9 1.22e-04 - 1.00e+00 1.00e+00h 1
234 1.9298310e-02 7.23e-05 4.77e-02 -7.3 6.01e-02 - 9.72e-01 1.00e+00f 1
235 1.8944265e-02 7.10e-05 2.58e-02 -7.3 3.37e-01 - 1.00e+00 4.94e-01h 1
236 1.8850512e-02 2.50e-04 3.20e-02 -7.3 1.22e+00 - 1.00e+00 3.01e-01h 1
237 1.8973382e-02 1.68e-04 1.62e-02 -5.6 2.88e-01 - 1.00e+00 9.69e-01f 1
238 1.8923222e-02 7.42e-04 1.98e-02 -5.7 6.98e-01 - 1.00e+00 8.97e-01h 1
239 1.8909911e-02 2.38e-06 1.07e-03 -5.7 5.19e-01 - 1.00e+00 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
240 1.8768536e-02 4.24e-04 5.99e-03 -6.4 4.13e-01 - 1.00e+00 1.00e+00h 1
241 1.8766483e-02 2.57e-07 1.07e-03 -7.2 7.18e-04 - 1.00e+00 1.00e+00h 1
242 1.8765610e-02 5.09e-09 1.59e-04 -9.3 1.60e-04 - 1.00e+00 1.00e+00h 1
243 1.8756402e-02 1.57e-07 3.71e-04 -7.3 5.44e-01 - 1.00e+00 3.20e-03h 1
244 1.8754641e-02 2.88e-07 4.23e-04 -6.9 4.17e-01 - 1.00e+00 1.27e-02h 1
245 1.8741111e-02 3.95e-05 2.52e-03 -7.8 2.87e-01 - 1.00e+00 2.34e-01h 1
246 1.8761846e-02 5.21e-04 2.20e-02 -6.3 1.57e-01 - 1.00e+00 1.00e+00h 1
247 1.8769273e-02 1.95e-07 3.84e-03 -6.3 9.67e-04 - 1.00e+00 1.00e+00h 1
248 1.8769143e-02 4.52e-10 9.17e-05 -6.3 9.16e-05 - 1.00e+00 1.00e+00h 1
249 1.8750172e-02 6.63e-06 3.61e-04 -7.1 8.99e-02 - 9.98e-01 3.30e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
250 1.8726350e-02 2.20e-05 1.46e-03 -7.0 6.38e-02 - 1.00e+00 7.63e-01h 1
251 1.8738789e-02 1.10e-05 9.57e-04 -6.5 1.18e-01 - 1.00e+00 1.00e+00h 1
252 1.8729265e-02 1.22e-05 3.00e-03 -6.6 5.94e-01 - 1.00e+00 7.57e-01h 1
253 1.8712596e-02 3.98e-05 4.25e-03 -7.3 5.69e-01 - 1.00e+00 5.28e-01h 1
254 1.8699983e-02 2.33e-05 4.81e-04 -7.2 6.96e-02 - 1.00e+00 9.97e-01h 1
255 1.8699721e-02 9.63e-10 8.13e-05 -9.3 1.04e-04 - 1.00e+00 1.00e+00h 1
256 1.8697901e-02 1.75e-07 5.53e-04 -10.3 1.94e-01 - 1.00e+00 2.53e-02h 1
257 1.8696388e-02 4.97e-06 3.16e-03 -9.8 2.32e+00 - 1.00e+00 1.39e-02h 1
258 1.8678973e-02 1.48e-03 3.59e-02 -7.5 1.62e+00 - 1.00e+00 5.64e-01h 1
259 1.8693706e-02 1.01e-03 1.99e-02 -6.8 1.48e+00 - 1.00e+00 6.67e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
260 1.8696262e-02 3.05e-08 2.48e-04 -6.5 1.68e-03 - 1.00e+00 1.00e+00h 1
261 1.8695635e-02 2.72e-09 9.45e-05 -8.5 9.48e-05 - 1.00e+00 1.00e+00h 1
262 1.8690779e-02 3.29e-07 9.21e-04 -8.3 1.03e-01 - 1.00e+00 5.43e-02h 1
263 1.8688481e-02 1.83e-05 9.47e-03 -9.0 9.50e-01 - 7.53e-01 1.14e-01h 1
264 1.8676480e-02 4.02e-03 9.55e-02 -6.7 1.18e+01 - 3.24e-01 1.68e-01h 1
265 1.8487857e-02 3.97e-03 2.74e-03 -7.0 9.58e-01 - 1.00e+00 1.00e+00h 1
266 1.8732211e-02 9.92e-05 8.11e-02 -6.2 2.09e-01 - 1.00e+00 1.00e+00h 1
267 1.8721073e-02 5.12e-05 5.33e-02 -6.3 6.84e-01 - 8.87e-01 6.56e-01h 1
268 1.8711300e-02 2.45e-04 6.54e-02 -6.3 2.25e+00 - 1.92e-01 1.06e-01f 1
269 1.8703295e-02 5.15e-04 8.01e-02 -6.3 6.35e+00 - 2.54e-01 5.77e-02f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
270 1.8695060e-02 4.42e-04 7.04e-02 -6.3 5.49e+00 - 2.62e-01 1.43e-01f 1
271 1.8694911e-02 3.98e-07 1.75e-01 -6.3 1.76e-01 - 1.00e+00 1.00e+00f 1
272 1.8694856e-02 1.32e-11 7.11e-05 -6.3 7.07e-05 - 1.00e+00 1.00e+00h 1
273 1.8674166e-02 8.23e-06 4.45e-03 -7.5 1.15e-01 - 1.00e+00 3.81e-01f 1
274 1.8643328e-02 2.86e-05 9.31e-03 -8.2 1.36e-01 - 1.00e+00 9.66e-01h 1
275 1.8645345e-02 6.35e-06 4.87e-03 -7.3 5.89e-02 - 1.00e+00 1.00e+00h 1
276 1.8644819e-02 2.24e-09 7.96e-05 -9.2 8.35e-05 - 1.00e+00 1.00e+00h 1
277 1.8643081e-02 1.80e-07 5.20e-04 -10.5 8.96e-02 - 1.00e+00 3.06e-02h 1
278 1.8642558e-02 7.43e-07 7.76e-04 -9.3 2.70e-01 - 1.00e+00 3.41e-02h 1
279 1.8631311e-02 1.08e-03 1.83e-02 -8.0 2.37e+00 - 1.00e+00 4.40e-01h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
280 1.8629545e-02 7.91e-04 9.53e-03 -7.7 3.94e-01 - 1.00e+00 2.68e-01h 1
281 1.8623293e-02 9.42e-05 1.32e-02 -7.7 5.08e-01 - 1.00e+00 1.00e+00h 1
282 1.8622005e-02 2.72e-08 1.05e-04 -8.4 1.58e-04 - 1.00e+00 1.00e+00h 1
283 1.8621728e-02 3.49e-08 8.55e-06 -9.5 1.45e-02 - 1.00e+00 1.91e-02h 1
284 1.8621424e-02 2.48e-07 6.97e-05 -8.9 8.88e-02 - 1.00e+00 6.40e-02h 1
285 1.8619949e-02 2.18e-05 3.67e-03 -7.6 7.10e-01 - 1.00e+00 1.35e-01f 1
286 1.8622162e-02 2.10e-04 1.20e-02 -6.9 7.53e-01 - 1.00e+00 4.28e-01h 1
287 1.8626788e-02 2.01e-06 1.05e-02 -7.0 1.26e-01 - 1.00e+00 1.00e+00h 1
288 1.8626827e-02 2.11e-11 1.82e-05 -7.0 1.86e-05 - 1.00e+00 1.00e+00h 1
289 1.8615731e-02 3.99e-06 3.26e-03 -8.4 6.03e-02 - 9.98e-01 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
290 1.8614460e-02 4.61e-06 2.46e-03 -7.9 6.49e-02 - 1.00e+00 9.96e-01h 1
291 1.8617335e-02 1.03e-06 1.69e-05 -7.4 1.62e-02 - 1.00e+00 1.00e+00h 1
292 1.8616759e-02 3.04e-09 7.90e-05 -9.4 8.01e-05 - 1.00e+00 1.00e+00h 1
293 1.8615250e-02 1.89e-08 1.54e-05 -10.7 5.98e-02 - 1.00e+00 3.09e-02h 1
294 1.8615000e-02 5.38e-08 5.34e-06 -10.3 4.26e-01 - 1.00e+00 2.22e-02h 1
295 1.8614756e-02 1.66e-07 2.67e-04 -10.3 3.92e-01 - 1.00e+00 3.82e-02h 1
296 1.8611397e-02 4.00e-05 5.66e-03 -9.1 3.04e-01 - 1.00e+00 9.05e-01h 1
297 1.8611053e-02 2.37e-09 3.53e-05 -9.6 7.73e-05 - 1.00e+00 1.00e+00h 1
298 1.8611046e-02 2.36e-09 5.30e-06 -11.0 7.27e-03 - 1.00e+00 9.15e-03h 1
299 1.8610128e-02 1.17e-05 3.74e-03 -10.4 1.90e+00 - 1.00e+00 5.84e-02h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
300 1.8607660e-02 1.44e-04 1.53e-02 -7.5 2.08e+00 - 1.00e+00 1.68e-01f 1
301 1.8614696e-02 1.80e-04 2.00e-02 -5.9 4.32e+00 - 1.00e+00 6.97e-02h 1
302 1.8618762e-02 1.25e-07 7.69e-04 -6.1 8.92e-04 - 1.00e+00 1.00e+00h 1
303 1.8626936e-02 5.63e-07 1.31e-03 -6.1 1.59e-02 - 1.00e+00 1.00e+00h 1
304 1.8682092e-02 2.02e-05 5.38e-04 -6.1 1.65e-01 - 1.00e+00 1.00e+00h 1
305 1.8682922e-02 2.17e-05 5.30e-03 -6.1 1.96e-01 - 1.00e+00 1.00e+00h 1
306 1.8681819e-02 1.77e-05 4.15e-03 -6.1 2.32e-01 - 1.00e+00 1.00e+00h 1
307 1.8679267e-02 9.71e-05 8.88e-03 -6.1 3.26e-01 - 1.00e+00 1.00e+00h 1
308 1.8678582e-02 9.76e-09 1.39e-05 -6.1 1.81e-04 - 1.00e+00 1.00e+00h 1
309 1.8677425e-02 8.89e-06 2.09e-03 -6.1 1.87e-02 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
310 1.8676884e-02 5.36e-06 1.29e-03 -6.1 2.14e-02 - 1.00e+00 1.00e+00h 1
311 1.8675435e-02 4.14e-05 4.56e-03 -6.1 4.11e-01 - 1.00e+00 3.69e-01h 2
312 1.8674092e-02 9.00e-05 8.20e-03 -6.1 4.94e-01 - 1.00e+00 1.00e+00h 1
313 1.8672740e-02 7.82e-05 9.11e-03 -6.1 4.89e-01 - 1.00e+00 3.62e-01h 2
314 1.8672312e-02 3.56e-09 1.32e-05 -6.1 1.54e-04 - 1.00e+00 1.00e+00h 1
315 1.8671810e-02 3.17e-06 1.46e-03 -6.1 1.14e-02 - 1.00e+00 1.00e+00h 1
316 1.8671585e-02 2.76e-06 9.61e-04 -6.1 1.48e-02 - 1.00e+00 1.00e+00h 1
317 1.8671310e-02 2.25e-06 9.23e-04 -6.1 1.66e-02 - 1.00e+00 1.00e+00h 1
318 1.8671305e-02 2.35e-13 4.37e-06 -6.1 4.70e-06 - 1.00e+00 1.00e+00h 1
319 1.8592433e-02 5.54e-06 6.04e-03 -8.2 1.17e-01 - 9.48e-01 9.98e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
320 1.8601276e-02 6.27e-07 4.30e-04 -6.9 1.11e-01 - 1.00e+00 9.59e-01h 1
321 1.8597617e-02 3.03e-06 6.31e-04 -7.3 6.76e-01 - 1.00e+00 3.50e-01h 1
322 1.8594246e-02 2.98e-06 8.04e-04 -11.0 6.36e-02 - 4.83e-01 3.02e-01h 1
323 1.8590137e-02 2.19e-05 4.22e-04 -8.0 2.16e-01 - 1.00e+00 5.87e-01h 1
324 1.8611815e-02 7.03e-05 9.70e-03 -6.6 4.31e-01 - 1.00e+00 1.00e+00h 1
325 1.8609932e-02 7.74e-05 1.12e-02 -6.7 9.12e-01 - 1.00e+00 2.32e-01h 2
326 1.8607142e-02 7.45e-05 1.04e-02 -6.7 5.20e-01 - 1.00e+00 4.06e-01h 1
327 1.8606919e-02 2.56e-08 5.32e-02 -6.7 5.32e-02 - 1.00e+00 1.00e+00f 1
328 1.8621619e-02 2.09e-06 8.08e-03 -5.3 2.46e-03 - 9.91e-01 1.00e+00h 1
329 1.8641377e-02 9.75e-07 2.65e-02 -5.4 5.12e-02 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
330 1.8978010e-02 1.07e-03 2.32e-02 -5.4 4.02e-01 - 1.00e+00 1.00e+00h 1
331 1.8970142e-02 1.46e-03 4.08e-02 -5.4 4.50e-01 - 1.00e+00 1.00e+00h 1
332 1.8981852e-02 1.45e-05 9.02e-03 -5.4 5.65e-02 - 1.00e+00 1.00e+00h 1
333 1.8983119e-02 9.01e-06 9.31e-03 -5.4 4.88e-02 - 1.00e+00 1.00e+00h 1
334 1.8983370e-02 3.42e-06 2.90e-03 -5.4 4.24e-02 - 1.00e+00 1.00e+00h 1
335 1.8983621e-02 1.53e-05 2.10e-03 -5.4 4.07e-02 - 1.00e+00 1.00e+00h 1
336 1.8983764e-02 2.28e-05 2.15e-03 -5.4 5.76e-02 - 1.00e+00 1.00e+00h 1
337 1.8983671e-02 2.98e-05 2.54e-03 -5.4 7.90e-02 - 1.00e+00 1.00e+00h 1
338 1.8983677e-02 1.65e-10 5.50e-06 -5.4 7.43e-05 - 1.00e+00 1.00e+00h 1
339 1.8714235e-02 5.78e-05 1.73e-02 -8.1 5.65e-02 - 9.73e-01 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
340 1.8632417e-02 3.16e-05 7.86e-03 -8.1 2.79e-01 - 1.00e+00 5.04e-01h 1
341 1.8616655e-02 2.68e-05 8.30e-03 -8.1 1.31e+00 - 1.00e+00 1.65e-01h 1
342 1.8612575e-02 3.35e-05 1.18e-02 -8.1 1.22e+01 - 9.54e-02 2.53e-02h 1
343 1.8609554e-02 4.73e-05 1.63e-02 -8.1 1.47e+01 - 1.12e-01 3.21e-02h 1
344 1.8602794e-02 1.25e-04 2.70e-02 -8.1 9.62e-01 - 1.82e-01 2.31e-01h 1
345 1.8601101e-02 1.31e-04 3.15e-02 -8.1 1.63e+00 - 1.00e+00 6.20e-02h 1
346 1.8602146e-02 8.08e-09 5.55e-03 -8.1 5.55e-03 - 1.00e+00 1.00e+00h 1
347 1.8596571e-02 5.62e-06 1.99e-03 -8.0 2.33e+00 - 9.96e-01 1.30e-02h 1
348 1.8585443e-02 4.22e-04 2.07e-02 -7.5 5.32e-01 - 1.00e+00 8.41e-01h 1
349 1.8579657e-02 2.88e-04 4.27e-04 -7.7 3.11e-01 - 1.00e+00 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
350 1.8583918e-02 7.88e-06 2.18e-04 -7.9 6.42e-02 - 1.00e+00 1.00e+00h 1
351 1.8582773e-02 1.02e-05 5.27e-04 -8.6 6.64e-02 - 1.00e+00 9.99e-01h 1
352 1.8582579e-02 3.72e-07 3.41e-04 -10.0 1.44e-02 - 1.00e+00 1.00e+00h 1
353 1.8582564e-02 1.88e-08 4.11e-07 -11.0 2.60e-02 - 1.00e+00 1.00e+00h 1
354 1.8582182e-02 8.56e-06 3.14e-04 -11.0 1.10e+00 - 1.00e+00 8.33e-01h 1
355 1.8581772e-02 8.58e-06 1.49e-04 -11.0 3.40e-01 - 1.00e+00 1.00e+00h 1
356 1.8585033e-02 2.32e-05 7.14e-04 -7.4 3.64e-01 - 1.00e+00 9.35e-01f 1
357 1.8583639e-02 2.58e-05 9.00e-04 -7.5 1.08e-01 - 1.00e+00 1.00e+00h 1
358 1.8583337e-02 1.99e-05 1.97e-03 -7.5 1.78e-01 - 1.00e+00 5.02e-01H 1
359 1.8582701e-02 6.90e-05 1.33e-03 -7.5 1.89e-01 - 1.00e+00 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
360 1.8582936e-02 1.74e-08 4.81e-07 -7.5 3.78e-03 - 1.00e+00 1.00e+00h 1
361 1.8579897e-02 1.88e-08 2.72e-07 -9.6 8.05e-03 - 1.00e+00 9.90e-01h 1
362 1.8579840e-02 6.06e-09 7.56e-08 -10.2 1.43e-02 - 1.00e+00 9.32e-01h 1
363 1.8579830e-02 4.38e-09 7.97e-08 -11.0 1.52e-02 - 1.00e+00 1.00e+00h 1
364 1.8579829e-02 4.41e-09 8.38e-08 -11.0 1.60e-02 - 1.00e+00 1.00e+00h 1
365 1.8579829e-02 1.11e-16 8.38e-08 -11.0 8.38e-08 - 1.00e+00 1.00e+00h 1
366 1.8579829e-02 1.11e-16 8.38e-08 -11.0 8.38e-08 - 1.00e+00 1.00e+00h 1
367 1.8579828e-02 5.12e-08 1.01e-07 -11.0 1.17e-02 - 1.00e+00 1.00e+00h 1
368 1.8579818e-02 7.22e-07 7.78e-04 -11.0 1.12e-01 - 1.00e+00 1.00e+00h 1
369 1.8579794e-02 1.53e-06 1.23e-03 -11.0 1.62e-01 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
370 1.8579683e-02 1.23e-05 2.45e-03 -11.0 3.54e-01 - 1.00e+00 1.00e+00h 1
371 1.8579597e-02 8.74e-11 4.60e-05 -9.6 5.00e-05 - 1.00e+00 1.00e+00h 1
372 1.8579595e-02 1.19e-10 7.92e-07 -11.0 5.87e-03 - 1.00e+00 1.70e-02h 1
373 1.8579064e-02 1.53e-04 1.01e-02 -11.0 1.77e+00 - 1.00e+00 1.74e-01h 1
374 1.8579610e-02 1.27e-05 1.96e-02 -9.9 4.17e-01 - 1.00e+00 1.00e+00H 1
375 1.8580271e-02 2.03e-06 1.47e-03 -8.2 4.20e-02 - 1.00e+00 1.00e+00h 1
376 1.8580271e-02 7.61e-13 1.54e-05 -8.2 1.55e-05 - 1.00e+00 1.00e+00h 1
377 1.8580148e-02 4.26e-07 7.78e-04 -8.2 9.73e-02 - 1.00e+00 1.00e+00h 1
378 1.8579846e-02 1.68e-06 1.01e-03 -8.2 1.29e+00 - 1.00e+00 2.15e-01h 1
379 1.8579735e-02 3.06e-06 1.20e-03 -8.2 5.21e-02 - 1.00e+00 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
380 1.8579452e-02 1.71e-05 3.10e-03 -8.2 1.23e-01 - 1.00e+00 6.50e-01h 1
381 1.8579035e-02 1.12e-06 5.34e-03 -8.7 2.36e-02 - 1.00e+00 1.00e+00h 1
382 1.8578852e-02 1.83e-06 7.00e-03 -8.9 2.89e-02 - 1.00e+00 1.00e+00h 1
383 1.8578647e-02 2.52e-06 8.57e-03 -9.9 3.52e-02 - 1.00e+00 1.00e+00h 1
384 1.8578630e-02 1.11e-12 7.25e-06 -11.0 7.58e-06 - 1.00e+00 1.00e+00h 1
385 1.8578627e-02 1.59e-09 1.80e-06 -11.0 2.60e-02 - 1.00e+00 3.26e-02h 1
386 1.8578550e-02 1.39e-06 1.53e-03 -11.0 2.66e-02 - 1.00e+00 1.00e+00h 1
387 1.8578473e-02 1.38e-06 1.54e-03 -11.0 2.67e-02 - 1.00e+00 1.00e+00h 1
388 1.8578473e-02 1.51e-13 1.28e-06 -11.0 3.29e-06 - 1.00e+00 1.00e+00h 1
389 1.8578473e-02 4.16e-14 1.80e-06 -11.0 1.80e-06 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
390 1.8577841e-02 1.01e-04 5.16e-04 -11.0 4.14e-01 - 1.00e+00 5.38e-01h 1
391 1.8577868e-02 3.73e-05 2.10e-04 -10.7 1.31e-01 - 1.00e+00 1.00e+00h 1
392 1.8578202e-02 3.00e-06 4.59e-04 -11.0 3.77e-02 - 1.00e+00 1.00e+00h 1
393 1.8578215e-02 4.64e-07 9.87e-05 -11.0 1.48e-02 - 1.00e+00 1.00e+00h 1
394 1.8578218e-02 1.40e-08 1.15e-08 -11.0 2.56e-03 - 1.00e+00 1.00e+00h 1
395 1.8578218e-02 2.28e-11 4.69e-10 -11.0 1.03e-04 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 395
(scaled) (unscaled)
Objective...............: 1.8578218284362957e-02 1.8578218284362957e-02
Dual infeasibility......: 4.6934814084267549e-10 4.6934814084267549e-10
Constraint violation....: 2.2769008900525023e-11 2.2769008900525023e-11
Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 1.0000000007918068e-11 1.0000000007918068e-11
Overall NLP error.......: 4.6934814084267549e-10 4.6934814084267549e-10
Number of objective function evaluations = 424
Number of objective gradient evaluations = 396
Number of equality constraint evaluations = 424
Number of inequality constraint evaluations = 424
Number of equality constraint Jacobian evaluations = 396
Number of inequality constraint Jacobian evaluations = 396
Number of Lagrangian Hessian evaluations = 0
Total seconds in IPOPT = 6.059
EXIT: Optimal Solution Found.
(Local) solution found
objective value = 0.018578218284362957
[7]:
# Value of the λ parameter, namely the constant control applied in the negative region
direct_sol.x3[1]
# Retrieving the corresponding switching time for plotting
N = length(direct_sol.t)
switch_index = findall(direct_sol.x1[2:N] .> no_control_bound)[1]; # switch from u = λ to u = 0
[8]:
# Plots
t = direct_sol.t
x1 = direct_sol.x1
x2 = direct_sol.x2
u = [direct_sol.x3[1:switch_index]; direct_sol.u[switch_index+1:end]]
p1 = direct_sol.p1
p2 = direct_sol.p2
x1_plot = plot(t, x1, xlabel = "t", ylabel = "x1", legend = false)
x2_plot = plot(t, x2, xlabel = "t", ylabel = "x2", legend = false)
p1_plot = plot(t, p1, xlabel = "t", ylabel = "p1", legend = false)
p2_plot = plot(t, p2, xlabel = "t", ylabel = "p2", legend = false)
u_plot = plot(t, u, xlabel = "t", ylabel = "u", legend = false, size=(800,400)) #, linetype=:steppre)
x1x2_plot = plot(x1, x2, xlabel = "x1", ylabel = "x2", legend = false)
plot(x1_plot, x2_plot, x1x2_plot, p1_plot, p2_plot, u_plot, layout = (3,2), size=(800,1200), left_margin=10mm)
[8]:
We can clearly identify four phases for this problem:
Constant control phase;
Zero control phase;
Bang-bang control, with successively u = -1 then u = +1.
Each of these phases have their own properties and can be identified mathematically, provided a first study from the direct method results. Exploiting this knowledge, we can state the different phases from the indirect solving point of view, since we do not provide the complete dynamics of the system.
In particular, for our system, the Hamiltonian for a given phase is $H = H_0 + H_1:nbsphinx-math:`cdot `u $; hence, the phases are identified this way:
Constant control to zero control: when x1 > -0.65;
Zero control to bang-bang (lower bound): x1’s sign changes;
Bang-bang (lower bound) to bang-bang (upper bound): pseudo-Hamiltonian’s sign changes.
Moreover, we also use the results from direct method as initial guess for the indirect solving, exploiting the same properties that we identified for the indirect solving.
Indirect solving¶
[9]:
# Imports for indirect method
using NLsolve
include("flow.jl");
[10]:
# Dynamics
function F0(x)
return [ x[2], 0.0 ]
end
function F1(x)
return [ 0., 1. ]
end
# Hamiltonians: permanent region
H0(x, p) = p' * F0(x)
H1(x, p) = p' * F1(x)
Hc(x, p, u) = H0(x, p) + u*H1(x,p) # pseudo-Hamiltonian
up(x, p) = M
um(x, p) = -M
Hp(x, p) = Hc(x, p, up(x, p))
Hm(x, p) = Hc(x, p, um(x, p))
# Hamiltonian: no control region
Hc0(x, p) = H0(x, p)
# Hamiltonians: control loss region
Hnc(x, p, λ) = H0(x, p) + λ*H1(x,p) # pseudo-Hamiltonian
Hnc(X, P) = Hnc(X[1:2], P[1:2], X[3])
# Flows
fp = Flow(Hp)
fm = Flow(Hm)
fc0 = Flow(Hc0)
fnc = Flow(Hnc);
[11]:
N = length(t); # number of time steps
x = [ [ x1[i], x2[i] ] for i in 1:N ];
p = -[ [ p1[i], p2[i] ] for i in 1:N ];
H1_span = H1.(x, p);
[12]:
H1_plot = plot(t, H1_span, xlabel = "t", ylabel = "H1", legend = false)
display(plot(u_plot, H1_plot, layout = (1,2), size=(800,400)))
[13]:
# Shooting function
# We know from direct method that the structure of the control is
# a constant region (λ) followed by bang-bang control (lower then upper bound)
function shootBλB0BmBp(p0, λ, t1, t2, t3, T, x0) # Bλ B0 B- B+ structure
# Phase 1: constant control
pλ0 = 0. # λ0 is free
X1, P1 = fnc(t0, [x0; λ], [p0; pλ0], t1) # augmented flow
# Phase 2: no control
jump_magnitude_1 = P1[2]*(-λ)/X1[2]
x2, p2 = fc0(t1, X1[1:2], [P1[1]-jump_magnitude_1, P1[2]], t2)
#x2, p2 = fc0(t1, X1[1:2], [P1[1], P1[2]], t2)
# Phase 3: bang-bang control, lower bound
jump_magnitude_2 = p2[2]*(-M)/x2[2]
x3, p3 = fm(t2, x2, [p2[1]-jump_magnitude_2, p2[2]], t3)
# Phase 4: bang-bang control, upper bound
xT, pT = fp(t3, x3, p3, T)
s = zeros(eltype(p0), 7)
s[1:2] = xT - [ x1f, x2f ] # target
s[3] = X1[1] - no_control_bound # switching to u = 0
s[4] = H1(x3, p3) # switching to u = +1
s[5] = Hp(xT, pT) - 1. # free final time
s[6] = x2[1] # crossing x1 = 0; u = -1
s[7] = P1[3] # ∫ ∂H∂λ = 0 pris sur [t0, t1]
return s
end;
[14]:
switch_index_1 = findall(x1[2:N] .> no_control_bound)[1]; # switch for u = 0
switch_index_2 = switch_index_1 + findall(x1[switch_index_1+1:N].*x1[switch_index_1:N-1].<0)[1]; # switch for u = -1
switch_index_3 = switch_index_2 + findall(H1_span[switch_index_2+1:N].*H1_span[switch_index_2:N-1].<0.)[1]; # switch for u = +1
ts1 = t[switch_index_1];
ts2 = t[switch_index_2];
ts3 = t[switch_index_3];
println("First switch time: ", ts1);
println("Second switch time: ", ts2);
println("Third switch time: ", ts3);
First switch time: 0.9289109142181479
Second switch time: 1.82066539186757
Third switch time: 3.102562453488614
[15]:
# Solve
x0 = [-1.; 0.]
S(ξ) = shootBλB0BmBp(ξ[1:2], ξ[3], ξ[4], ξ[5], ξ[6], ξ[7], x0)
p0 = [p1[1], p2[1]]
tf = direct_sol.t[end]
λs = direct_sol.x3[1]
# ξ contains shoot function parameters
ξ_guess = [p0..., λs, ts1, ts2, ts3, tf] # initial guess (from direct method)
println("Initial value of shooting:\n", S(ξ_guess), "\n\n")
indirect_sol = nlsolve(S, ξ_guess; xtol=1e-8, method=:trust_region, show_trace=true)
println(indirect_sol)
println("Final value of shooting:\n", S(indirect_sol.zero), "\n\n")
# Retrieves solution
if indirect_sol.f_converged || indirect_sol.x_converged
p0 = indirect_sol.zero[1:2]
λ = indirect_sol.zero[3]
t1 = indirect_sol.zero[4]
t2 = indirect_sol.zero[5]
t3 = indirect_sol.zero[6]
T = indirect_sol.zero[7]
else
error("Not converged")
end;
Initial value of shooting:
[-0.028780650249423726, 0.06655286031403253, -0.008453985681637022, -0.00036371023282586406, -1.0057241217301314, -0.0026856381903773052, 7.515869035458086e-5]
Iter f(x) inf-norm Step 2-norm
------ -------------- --------------
0 1.005724e+00 NaN
1 9.853084e+00 4.578126e-16
2 4.959772e+00 0.000000e+00
3 8.294626e-01 2.220446e-16
4 4.190181e-01 5.264489e-01
5 1.793475e+00 6.206335e-17
6 4.720322e-01 5.611772e-01
7 3.493486e-01 2.257742e-01
8 2.958390e-01 4.589893e-01
9 2.841551e-01 4.749293e-01
10 4.458390e-01 0.000000e+00
11 2.392055e-01 6.818696e-01
12 2.907312e-01 2.482534e-16
13 1.406382e-01 6.969745e-01
14 8.160577e-03 6.522708e-01
15 9.046162e-05 1.030036e-02
16 2.992035e-09 3.524767e-05
Results of Nonlinear Solver Algorithm
* Algorithm: Trust-region with dogleg and autoscaling
* Starting Point: [-0.013415443489824715, -0.006311786485571456, 0.791646117292151, 0.9289109142181479, 1.82066539186757, 3.102562453488614, 3.7156436568725915]
* Zero: [2.6745514341370815, 1.2518144407624485, 0.7988404408201301, 0.9360930021555836, 1.8053222184490503, 3.081877799801171, 3.6106444346058764]
* Inf-norm of residuals: 0.000000
* Iterations: 16
* Convergence: true
* |x - x'| < 1.0e-08: false
* |f(x)| < 1.0e-08: true
* Function Calls (f): 17
* Jacobian Calls (df/dx): 11
Final value of shooting:
[-3.548628719089967e-11, -5.681065585522522e-11, -5.042588568926476e-11, 1.654177993106719e-9, 2.992035286197847e-9, -1.3948028836773354e-10, -5.803724346470672e-11]
[16]:
# plots
pλ0 = 0.
ode_sol = fnc((t0, t1), [x0; λ], [p0; pλ0], saveat=0.1)
tt1 = ode_sol.t
xx1 = [ ode_sol[1:2, j] for j in 1:size(tt1, 1) ]
pp1 = [ ode_sol[4:5, j] for j in 1:size(tt1, 1) ]
uu1 = λ.*ones(length(tt1))
ν_1 = ode_sol[5, end]*(-λ)/ode_sol[2, end] # jump
ode_sol = fc0((t1, t2), xx1[end], pp1[end]+[ν_1, 0], saveat=0.1)
tt2 = ode_sol.t
xx2 = [ ode_sol[1:2, j] for j in 1:size(tt2, 1) ]
pp2 = [ ode_sol[3:4, j] for j in 1:size(tt2, 1) ]
uu2 = zeros(size(tt2, 1))
ν_2 = ode_sol[4, end]*(-M)/ode_sol[2, end] # jump
ode_sol = fm((t2, t3), xx2[end], pp2[end]+[ν_2, 0], saveat=0.1)
tt3 = ode_sol.t
xx3 = [ ode_sol[1:2, j] for j in 1:size(tt3, 1) ]
pp3 = [ ode_sol[3:4, j] for j in 1:size(tt3, 1) ]
uu3 = um.(xx3, pp3)
ode_sol = fp((t3, T), xx3[end], pp3[end], saveat=0.1)
tt4 = ode_sol.t
xx4 = [ ode_sol[1:2, j] for j in 1:size(tt4, 1) ]
pp4 = [ ode_sol[3:4, j] for j in 1:size(tt4, 1) ]
uu4 = up.(xx4, pp4)
t = [ tt1 ; tt2 ; tt3 ; tt4 ]
x = [ xx1 ; xx2 ; xx3 ; xx4 ]
p = [ pp1 ; pp2 ; pp3 ; pp4 ]
u = [ uu1 ; uu2 ; uu3 ; uu4 ]
m = length(t)
x1 = [ x[i][1] for i=1:m ]
x2 = [ x[i][2] for i=1:m ]
p1 = [ p[i][1] for i=1:m ]
p2 = [ p[i][2] for i=1:m ];
x1_plot = plot(t, x1, xlabel = "t", ylabel = "x1", legend = false)
x2_plot = plot(t, x2, xlabel = "t", ylabel = "x2", legend = false)
p1_plot = plot(t, p1, xlabel = "t", ylabel = "p1", legend = false)
p2_plot = plot(t, p2, xlabel = "t", ylabel = "p2", legend = false)
u_plot = plot(t, u, xlabel = "t", ylabel = "u", legend = false, size=(800,400)) #, linetype=:steppre)
x1x2_plot = plot(x1, x2, xlabel = "x1", ylabel = "x2", legend = false)
plot(x1_plot, x2_plot, x1x2_plot, u_plot, p1_plot, p2_plot, layout = (3,2), size=(800,1200), left_margin=10mm)
[16]: