Notebook source code: examples/control-loss/double-integrator-with-control-loss-and-no-control.ipynb
Run the notebook yourself on binder Binder badge

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:

\[\begin{split}\left\{ \begin{array}{l} \displaystyle \min T, \\[0.5em] \dot{x}_1(t) = x_2(t) \\ \dot{x}_2(t) =\mathbb{1}_{\mathbb{R}⁻\setminus[-0.65, 0]}(x_1(t))\lambda + (1 - \mathbb{1}_{\mathbb{R}⁻}x_1(t)) u(t), \\[0.5em] \lambda,u(t) \in [-1, 1], \\[0.5em] x(0) = x_0,\quad x(T) = 0_{\mathrm{R}^2}. \end{array} \right.\end{split}\]

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.

Thumbnail

[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]:
../_images/control-loss_double-integrator-with-control-loss-and-no-control_8_0.svg

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)))
../_images/control-loss_double-integrator-with-control-loss-and-no-control_14_0.svg
[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]:
../_images/control-loss_double-integrator-with-control-loss-and-no-control_18_0.svg