Note
Go to the end to download the full example code.
How scaling can help convergence#
This example showcases a optimal control task where the scaling of the NLP variables can really impact the performance of the solver. This is due to the fact that the variables have widely different scales, making it challenging for the solver to find a good numerical solutions.
We’ll skip over some of the formalities of the library, as we assume you are already a
bit familiar with it. If this is not the case, for more introductory examples on
csnlp, see the other Introductory examples.
The problem reproduces the example from this CasADi blog post, where an optimal control problem for a 1-dimensional rocket task is tackled. The goal is to attain a given altitude \(Y_f\) at time \(T\) with minimal fuel consumption. In mathematical terms, the rocket state at time \(t\) is given by \(x_t = [y_t, v_t, m_t]^\top\), i.e., the altitude, speed, and mass of the rocket. The continuous-time dynamics of the rocket are given by
\[\begin{split}f(x_t, u_t) = \begin{bmatrix}
v_t \\ \frac{u_t}{m_t} - g \\ -\alpha u_t
\end{bmatrix},\end{split}\]
where \(u\) is the thrust, and other symbols are constants. By Euler-discretizing the dynamics, we formulate the final optimal control problem
\[\begin{split}\begin{aligned}
\min_{\substack{u_0,\dots,u_{N-1} \\ x_0,\dots,x_N}} \quad & m_0 - m_N \\
\text{s.t.} \quad & y_N = Y_f \\
& x_0 = [0, 0, m_0]^\top \\
& x_{t+1} = x_t + f(x_t, u_t) \Delta t, \quad t = 0, \ldots, N-1.
\end{aligned}\end{split}\]
Without scaling#
First, we’ll address the problem without scaling. As usual, we’ll import the necessary
classes and modules from csnlp.
import casadi as cs
import matplotlib.pyplot as plt
import numpy as np
from csnlp import Nlp, scaling, wrappers
Let us define constants of the problem, the dynamics of the rocket, as well as the options for the IPOPT solver.
N = 100 # number of control intervals
T = 100 # Time horizon [s]
K = 3 # Nlp multistarts
dt = T / N
m0 = 500000 # start mass [kg]
yT = 100000 # final height [m]
g = 9.81 # gravity 9.81 [m/s^2]
alpha = 1 / (300 * g) # kg/(N*s)
seed = 69
opts = {"print_time": False, "ipopt": {"sb": "yes", "print_level": 5}}
def get_dynamics() -> cs.Function:
x, u = cs.SX.sym("x", 3), cs.SX.sym("u")
x_next = x + cs.vertcat(x[1], u / x[2] - g, -alpha * u) * dt
return cs.Function("F", [x, u], [x_next], ["x", "u"], ["x+"])
F = get_dynamics()
Now, we can build the optimization problem and solve it.
nlp = Nlp[cs.SX]()
x, _, _ = nlp.variable("x", (3, N))
u, _, _ = nlp.variable("u", (1, N - 1), lb=0)
nlp.constraint("dynamics", x[:, 1:], "==", F(x[:, :-1], u))
nlp.constraint("initial_conditions", x[:, 0], "==", [0, 0, m0])
nlp.constraint("final_conditions", x[0, -1], "==", yT)
nlp.minimize(x[2, 0] - x[2, -1])
nlp.init_solver(opts)
x_guess = cs.repmat([0, 0, 1e5], 1, N)
sol = nlp(vals0={"x": x_guess})
This is Ipopt version 3.14.11, running with linear solver MUMPS 5.4.1.
Number of nonzeros in equality constraint Jacobian...: 994
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 198
Total number of variables............................: 399
variables with only lower bounds: 99
variables with lower and upper bounds: 0
variables with only upper bounds: 0
Total number of equality constraints.................: 301
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 0.0000000e+00 4.00e+05 1.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 7.6368935e+04 3.63e+01 8.25e+08 -1.0 4.54e+06 - 2.18e-09 1.00e+00h 1
2 7.9629768e+04 3.60e+01 8.16e+08 -1.0 1.39e+07 0.0 9.42e-02 1.00e-02h 1
3 9.7228569e+04 3.37e+01 7.55e+08 -1.0 1.20e+07 0.4 1.77e-01 6.33e-02h 1
4 2.4523820e+05 3.93e+01 1.13e+09 -1.0 1.34e+07 -0.1 2.24e-01 1.00e+00H 1
5 3.1182841e+05 1.95e+00 2.35e+08 -1.0 6.31e+06 -0.5 5.10e-01 1.00e+00h 1
6 3.0175491e+05 1.26e-01 1.05e+07 -1.0 5.42e+05 -1.0 9.78e-01 1.00e+00f 1
7 3.0162314e+05 1.38e-05 3.01e+03 -1.0 6.24e+03 -1.5 9.90e-01 1.00e+00f 1
8 3.0162312e+05 5.55e-11 6.42e-02 -1.0 5.88e+00 -2.0 9.94e-01 1.00e+00f 1
9 3.0162312e+05 5.09e-11 7.49e-04 -1.7 2.05e-01 -2.4 1.00e+00 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 3.0162311e+05 5.28e-11 3.33e-04 -3.8 2.73e-01 -2.9 1.00e+00 1.00e+00f 1
11 3.0162311e+05 5.05e-11 3.33e-04 -5.7 8.19e-01 -3.4 1.00e+00 1.00e+00f 1
12 3.0162309e+05 4.91e-11 3.36e-04 -5.7 2.48e+00 -3.9 1.00e+00 1.00e+00f 1
13 3.0162307e+05 7.46e-11 3.33e-04 -5.7 7.37e+00 -4.3 1.00e+00 4.02e-01f 1
14 3.0162293e+05 7.80e-11 3.33e-04 -5.7 2.21e+01 -4.8 1.00e+00 1.00e+00f 1
15 3.0162249e+05 7.02e-10 3.33e-04 -5.7 6.63e+01 -5.3 1.00e+00 1.00e+00f 1
16 3.0162118e+05 6.31e-09 3.33e-04 -5.7 1.99e+02 -5.8 1.00e+00 1.00e+00f 1
17 3.0161726e+05 5.68e-08 3.33e-04 -5.7 5.97e+02 -6.3 1.00e+00 1.00e+00f 1
18 3.0160550e+05 5.11e-07 3.33e-04 -5.7 1.79e+03 -6.7 1.00e+00 1.00e+00f 1
19 3.0157020e+05 4.60e-06 3.33e-04 -5.7 5.37e+03 -7.2 1.00e+00 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20 3.0146429e+05 4.14e-05 3.33e-04 -5.7 1.61e+04 -7.7 1.00e+00 1.00e+00f 1
21 3.0114643e+05 3.72e-04 3.33e-04 -5.7 4.84e+04 -8.2 1.00e+00 1.00e+00f 1
22 3.0019168e+05 3.32e-03 3.33e-04 -5.7 1.45e+05 -8.6 1.00e+00 1.00e+00f 1
23 2.9991873e+05 3.28e-03 3.33e-04 -5.7 4.35e+05 -9.1 1.00e+00 9.49e-02f 1
24 2.9991873e+05 3.28e-03 3.26e-04 -5.7 1.28e+06 -9.6 1.00e+00 7.96e-07f 1
25 2.9819434e+05 1.36e-02 3.26e-04 -5.7 3.85e+06 -10.1 1.00e+00 6.86e-02f 1
26 2.9819431e+05 1.36e-02 3.19e-04 -5.7 1.41e+06 -9.6 1.00e+00 3.67e-06f 1
27 2.9648437e+05 2.26e-02 3.19e-04 -5.7 4.26e+06 -10.1 1.00e+00 6.32e-02f 1
28 2.9647485e+05 2.26e-02 3.13e-04 -5.7 4.14e+07 -11.1 6.36e-01 3.02e-05f 1
29 2.9423555e+05 3.50e-02 3.05e-04 -5.7 4.80e+07 -11.1 1.00e+00 5.99e-03f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
30 2.9377723e+05 3.36e-02 3.01e-04 -5.7 1.90e+06 -9.8 1.00e+00 3.99e-02f 1
31 2.9212769e+05 3.89e-02 2.95e-04 -5.7 5.67e+06 -10.3 1.00e+00 4.72e-02f 1
32 2.8806498e+05 7.91e-02 2.85e-04 -5.7 6.62e+07 -11.2 3.64e-01 8.41e-03f 1
33 2.8595977e+05 8.84e-02 2.48e-02 -5.7 1.90e+07 -10.8 1.00e+00 1.84e-02f 1
34 2.8086039e+05 1.47e-01 1.06e+00 -5.7 7.37e+07 -11.3 2.93e-01 1.10e-02f 1
35 2.7834641e+05 1.56e-01 1.03e+00 -5.7 1.97e+07 -10.9 1.00e+00 2.37e-02f 1
36 2.7291782e+05 2.02e-01 1.02e+00 -5.7 8.37e+07 -11.3 3.07e-01 1.17e-02f 1
37 2.7037377e+05 1.98e-01 9.93e-01 -5.7 1.99e+07 -10.9 1.00e+00 2.77e-02f 1
38 2.6508639e+05 2.12e-01 9.80e-01 -5.7 8.63e+07 -11.4 2.79e-01 1.34e-02f 1
39 2.6283824e+05 2.06e-01 9.51e-01 -5.7 1.98e+07 -11.0 1.00e+00 2.94e-02f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
40 2.5834335e+05 1.94e-01 5.74e+00 -5.7 1.93e+07 -11.4 1.00e+00 6.08e-02f 1
41 2.5651414e+05 1.88e-01 5.83e+00 -5.7 1.84e+07 -11.0 1.00e+00 3.08e-02f 1
42 2.5273807e+05 1.85e-01 5.75e+00 -5.7 8.18e+07 -11.5 3.29e-01 1.48e-02f 1
43 2.5092631e+05 1.78e-01 5.64e+00 -5.7 1.89e+07 -11.1 1.00e+00 3.56e-02f 1
44 2.4840286e+05 1.69e-01 5.36e+00 -5.7 2.06e+07 -11.1 1.00e+00 4.91e-02f 1
45 2.4711990e+05 1.55e-01 4.92e+00 -5.7 6.74e+06 -10.7 1.00e+00 8.32e-02f 1
46 2.4520942e+05 1.49e-01 4.66e+00 -5.7 2.08e+07 -11.2 1.00e+00 4.28e-02f 1
47 2.4393966e+05 1.35e-01 4.24e+00 -5.7 7.01e+06 -10.7 1.00e+00 9.01e-02f 1
48 2.4242547e+05 1.29e-01 4.00e+00 -5.7 1.79e+07 -11.2 1.00e+00 4.42e-02f 1
49 2.4146529e+05 1.20e-01 3.71e+00 -5.7 7.34e+06 -10.8 1.00e+00 7.25e-02f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
50 2.3957406e+05 1.19e-01 3.67e+00 -5.7 1.01e+08 -11.7 3.02e-01 1.05e-02f 1
51 2.3863371e+05 1.16e-01 3.63e+00 -5.7 2.42e+07 -11.3 1.00e+00 2.49e-02f 1
52 2.3685682e+05 1.14e-01 3.62e+00 -5.7 9.07e+07 -11.8 3.56e-01 1.25e-02f 1
53 2.3596666e+05 1.11e-01 3.53e+00 -5.7 2.56e+07 -11.4 1.00e+00 2.48e-02f 1
54 2.3482710e+05 1.08e-01 3.42e+00 -5.7 2.77e+07 -11.4 1.00e+00 3.08e-02f 1
55 2.3324214e+05 1.07e-01 2.85e+00 -5.7 1.29e+08 -11.9 5.52e-01 9.24e-03f 1
56 2.3212784e+05 1.10e-01 2.75e+00 -5.7 2.71e+07 -11.5 1.00e+00 3.49e-02f 1
57 2.3049398e+05 1.26e-01 2.71e+00 -5.7 9.60e+07 -11.9 3.56e-01 1.50e-02f 1
58 2.2861467e+05 1.56e-01 2.70e+00 -5.7 1.30e+09 -12.4 3.31e-02 1.30e-03f 1
59 2.2733089e+05 1.70e-01 2.65e+00 -5.7 8.38e+07 -12.0 1.00e+00 1.64e-02f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
60 2.2546369e+05 2.13e-01 2.63e+00 -5.7 3.03e+08 -12.5 1.03e-01 6.96e-03f 1
61 2.2413823e+05 2.38e-01 2.62e+00 -5.7 3.18e+08 -12.5 2.43e-01 5.41e-03f 1
62 2.2247982e+05 2.89e-01 2.74e+00 -5.7 1.35e+09 -13.0 2.82e-02 1.99e-03f 1
63 2.2165986e+05 3.02e-01 2.75e+00 -5.7 2.43e+08 -12.6 1.00e+00 5.86e-03f 1
64 2.1974052e+05 4.04e-01 2.73e+00 -5.7 4.24e+08 -13.1 9.21e-02 1.01e-02f 1
65 2.1886192e+05 4.30e-01 2.71e+00 -5.7 6.67e+08 - 1.38e-01 3.80e-03f 1
66 2.1765979e+05 4.87e-01 2.69e+00 -5.7 4.42e+08 - 2.67e-01 9.11e-03f 1
67 2.1669837e+05 5.31e-01 2.66e+00 -5.7 2.99e+08 - 5.85e-01 1.29e-02f 1
68 2.1573320e+05 5.86e-01 2.60e+00 -5.7 2.23e+08 - 2.82e-01 2.08e-02f 1
69 2.1510232e+05 6.09e-01 2.55e+00 -5.7 1.79e+08 - 1.00e+00 1.98e-02f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
70 2.1416296e+05 6.80e-01 2.44e+00 -5.7 1.48e+08 - 4.10e-01 4.20e-02f 1
71 2.1370292e+05 6.90e-01 2.37e+00 -5.7 1.20e+08 - 1.00e+00 3.08e-02f 1
72 2.1313653e+05 7.12e-01 2.25e+00 -5.7 1.05e+08 - 6.61e-01 4.97e-02f 1
73 2.1277848e+05 7.13e-01 2.15e+00 -5.7 8.99e+07 - 1.00e+00 4.36e-02f 1
74 2.1249837e+05 7.05e-01 2.06e+00 -5.7 8.07e+07 - 9.58e-01 4.25e-02f 1
75 2.1217735e+05 6.96e-01 1.94e+00 -5.7 7.39e+07 - 1.00e+00 5.88e-02f 1
76 2.1177507e+05 6.92e-01 1.75e+00 -5.7 6.43e+07 - 1.00e+00 9.94e-02f 1
77 2.1173099e+05 6.82e-01 1.72e+00 -5.7 5.33e+07 - 1.00e+00 1.63e-02f 1
78 2.1138020e+05 6.72e-01 1.49e+00 -5.7 5.25e+07 - 1.00e+00 1.34e-01f 1
79 2.1136177e+05 6.65e-01 1.47e+00 -5.7 4.22e+07 - 1.00e+00 1.12e-02f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
80 2.1102724e+05 6.33e-01 1.17e+00 -5.7 4.18e+07 - 1.00e+00 2.07e-01f 1
81 2.1102586e+05 6.32e-01 1.17e+00 -5.7 3.11e+07 - 1.00e+00 1.58e-03f 1
82 2.1072807e+05 5.83e-01 7.65e-01 -5.7 3.11e+07 - 1.00e+00 3.44e-01f 1
83 2.1072804e+05 5.83e-01 7.65e-01 -5.7 1.96e+07 - 1.00e+00 6.89e-05f 1
84 2.1051661e+05 4.33e-01 2.63e-01 -5.7 1.96e+07 - 1.00e+00 6.56e-01f 1
85 2.1051661e+05 4.33e-01 2.63e-01 -5.7 7.21e+06 - 1.00e+00 6.65e-06h 1
86 2.1051530e+05 8.09e-02 5.86e-04 -5.7 6.88e+06 - 1.00e+00 9.98e-01h 1
87 2.1052797e+05 2.30e-05 1.13e-08 -5.7 1.45e+05 - 1.00e+00 1.00e+00h 1
88 2.1052797e+05 5.13e-11 1.84e-11 -5.7 1.33e+01 - 1.00e+00 1.00e+00h 1
89 2.1052797e+05 5.64e-11 2.51e-14 -8.6 9.89e-01 - 1.00e+00 1.00e+00f 1
Number of Iterations....: 89
(scaled) (unscaled)
Objective...............: 2.1052796979582077e+05 2.1052796979582077e+05
Dual infeasibility......: 2.5059110602548840e-14 2.5059110602548840e-14
Constraint violation....: 5.6388671509921551e-11 5.6388671509921551e-11
Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 2.5065590837648647e-09 2.5065590837648647e-09
Overall NLP error.......: 2.5065590837648647e-09 2.5065590837648647e-09
Number of objective function evaluations = 91
Number of objective gradient evaluations = 90
Number of equality constraint evaluations = 91
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 90
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 89
Total seconds in IPOPT = 0.302
EXIT: Optimal Solution Found.
We can visualize the optimal thrust strategy and the predicted evolution of the states
_, axs = plt.subplots(4, 1, sharex=True, constrained_layout=True)
time = np.linspace(0, T, N)
u_opt = sol.value(u).full().flatten()
x_opt = sol.value(x).full()
axs[0].plot(time[:-1], u_opt)
for x_, ax in zip(x_opt, axs[1:]):
ax.plot(time, x_)
axs[0].set_ylabel("Thrust [N]")
axs[1].set_ylabel("Height [m]")
axs[2].set_ylabel("Speed [m/s]")
axs[3].set_ylabel("Mass [kg]")
axs[3].set_xlabel("Time [s]")
plt.show()
Though it is much more interesting to see how the solver behaved during its iterations. In particular, we can plot the convergence of primal and dual feasibility, which is pretty bad. This is due to the fact that the altitude and mass of the rocket have larger scales than the speed, which makes the solver struggle to find a good solution.
iter_stats = sol.stats["iterations"]
_, ax = plt.subplots(constrained_layout=True)
ax.semilogy(iter_stats["inf_pr"], label="Primal")
ax.semilogy(iter_stats["inf_du"], label="Dual")
ax.legend()
ax.set_xlabel("Iteration")
ax.set_ylabel("Feasibility")
plt.show()
With (automatic) scaling#
How can we then improve convergence? One way is to scale the variables of the NLP.
csnlp provides the scaling wrapper csnlp.wrappers.NlpScaling that can
do just this.
First, let us define the scales of each state and build an instance of
csnlp.core.scaling.Scaler, which contains all the information needed for
scaling. We can do so by calling the csnlp.core.scaling.Scaler.register.
x_nominal = cs.DM([1e5, 2e3, 3e5])
u_nominal = 1e8
scaler = scaling.Scaler()
scaler.register("x", scale=x_nominal)
scaler.register("u", scale=u_nominal)
Then, we wrap the NLP isntance in the scaling wrapper and solve the problem again. Since this wrapper is non-retroactive, we need to build the NLP again from scratch. The solver should now converge much faster.
nlp = Nlp[cs.SX]()
nlp = wrappers.NlpScaling[cs.SX](nlp, scaler=scaler) # added this line!
x, _, _ = nlp.variable("x", (3, N))
u, _, _ = nlp.variable("u", (1, N - 1), lb=0)
nlp.constraint("dynamics", x[:, 1:], "==", F(x[:, :-1], u))
nlp.constraint("initial_conditions", x[:, 0], "==", [0, 0, m0])
nlp.constraint("final_conditions", x[0, -1], "==", yT)
nlp.minimize(x[2, 0] - x[2, -1])
nlp.init_solver(opts)
sol_scaled = nlp(vals0={"x": x_guess})
This is Ipopt version 3.14.11, running with linear solver MUMPS 5.4.1.
Number of nonzeros in equality constraint Jacobian...: 994
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 198
Total number of variables............................: 399
variables with only lower bounds: 99
variables with lower and upper bounds: 0
variables with only upper bounds: 0
Total number of equality constraints.................: 301
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 0.0000000e+00 4.00e+05 7.70e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 -8.6851002e+03 3.97e+05 7.64e+00 -1.0 1.88e+01 - 1.25e-02 7.64e-03h 1
2 3.2497635e+04 2.74e+05 2.31e+01 -1.0 1.32e+00 - 1.00e+00 3.09e-01h 1
3 1.7236698e+05 7.28e+01 1.51e+01 -1.0 9.15e-01 - 2.74e-01 1.00e+00h 1
4 2.3286737e+05 2.93e+00 1.07e+00 -1.0 2.24e-01 - 1.00e+00 1.00e+00h 1
5 2.3637015e+05 2.40e+00 9.08e-02 -1.0 2.79e-01 - 1.00e+00 1.00e+00f 1
6 2.1210625e+05 5.58e+00 4.29e-01 -2.5 5.11e-01 - 6.21e-01 1.00e+00f 1
7 2.1137822e+05 4.35e-01 2.82e-02 -2.5 6.52e-02 - 9.86e-01 1.00e+00h 1
8 2.1057760e+05 2.23e-01 1.57e-02 -3.8 3.95e-02 - 9.86e-01 1.00e+00h 1
9 2.1056890e+05 2.51e-02 1.77e-03 -3.8 1.79e-02 - 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.1056816e+05 1.97e-03 1.44e-04 -3.8 6.45e-03 - 1.00e+00 1.00e+00h 1
11 2.1052848e+05 8.16e-04 5.96e-05 -5.7 2.92e-03 - 1.00e+00 1.00e+00h 1
12 2.1052845e+05 2.61e-06 1.87e-07 -5.7 2.46e-04 - 1.00e+00 1.00e+00h 1
13 2.1052796e+05 1.02e-07 7.46e-09 -8.6 3.11e-05 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 13
(scaled) (unscaled)
Objective...............: 7.0175985088253512e+01 2.1052795526476053e+05
Dual infeasibility......: 7.4646067105277325e-09 2.2393820131583198e-05
Constraint violation....: 5.0998830936066500e-09 1.0199766187213299e-07
Variable bound violation: 9.7787537744512781e-09 9.7787537744512781e-09
Complementarity.........: 4.3175866741765892e-09 1.2952760022529768e-05
Overall NLP error.......: 7.4646067105277325e-09 2.2393820131583198e-05
Number of objective function evaluations = 14
Number of objective gradient evaluations = 14
Number of equality constraint evaluations = 14
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 14
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 13
Total seconds in IPOPT = 0.024
EXIT: Optimal Solution Found.
We can again plot the optimal values, both for the scaled (solid, blue) and unscaled variables (dashed, orange). As you can see, there is no difference (aside of course from the scale) between the two solutions.
_, axs = plt.subplots(4, 1, sharex=True, constrained_layout=True)
u_opt_scaled = sol_scaled.value(u).full().flatten()
x_opt_scaled = sol_scaled.value(x).full()
axs[0].plot(time[:-1], u_opt_scaled)
for x_, ax in zip(x_opt_scaled, axs[1:]):
ax.plot(time, x_)
axs_twin = [ax.twinx() for ax in axs]
u_opt_unscaled = sol_scaled.value(nlp.unscale(u)).full().flatten()
x_opt_unscaled = sol_scaled.value(nlp.unscale(x)).full()
axs_twin[0].plot(time[:-1], u_opt_unscaled, "--", color="C1")
for x_, ax in zip(x_opt_unscaled, axs_twin[1:]):
ax.plot(time, x_, "--", color="C1")
for ax in axs:
ax.spines["left"].set_color("C0")
ax.tick_params(axis="y", colors="C0")
for ax in axs_twin:
ax.spines["left"].set_visible(False)
ax.spines["right"].set_color("C1")
ax.tick_params(axis="y", colors="C1")
axs[0].set_ylabel("Thrust [N]")
axs[1].set_ylabel("Height [m]")
axs[2].set_ylabel("Speed [m/s]")
axs[3].set_ylabel("Mass [kg]")
axs[3].set_xlabel("Time [s]")
plt.show()
Finally, we can plot the convergence of the solver. For comparison, we plot again the convergence of the unscaled problem (solid, faded), and on top of those also the ones of the scaled problem (dashed). As you can see, the scaled problem converges much faster (i.e., with less iterations and falling steeper).
iter_stats_scaled = sol_scaled.stats["iterations"]
_, ax = plt.subplots(constrained_layout=True)
ax.semilogy(iter_stats["inf_pr"], alpha=0.2, label="Primal", color="C0")
ax.semilogy(iter_stats["inf_du"], alpha=0.2, label="Dual", color="C1")
ax.semilogy(iter_stats_scaled["inf_pr"], "--", lw=2, color="C0")
ax.semilogy(iter_stats_scaled["inf_du"], "--", lw=2, color="C1")
ax.legend()
ax.set_xlabel("Iteration")
ax.set_ylabel("Feasibility")
plt.show()
Total running time of the script: (0 minutes 1.290 seconds)