Pyomo.DAE Example: Race Car - Optimization for Decision Science

Adapted from Pyomo.DAE: Racing Example Revisited and examples/dae/car_example.py

Learning Objectives¶

Introduces syntax for Pyomo.DAE
Shows time-scaling modeling trick
Practice using the units feature in Pyomo, and learn where in the workflow to check them (see Units in Pyomo.dae Models)

Install Packages and Load Modules¶

# Imports
import sys

if "google.colab" in sys.modules:
    !wget "https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/helper.py"
    import helper

    helper.easy_install()
else:
    sys.path.insert(0, "../")
    import helper
helper.set_plotting_style()

import pyomo.environ as pyo
import matplotlib.pyplot as plt
import pyomo.dae as dae
from pyomo.environ import units
from pyomo.util.check_units import assert_units_consistent

Problem Statement and Optimal Control Formulation¶

You are a race car driver with a simple goal. Drive distance $L$ in the minimal amount of time but come to a complete stop at the finish line.

Mathematically, you want to solve the following optimal control problem:

\begin{align*} \min_{u,t_f} \quad & t_f \\ \mathrm{s.t.} \quad & \frac{dx}{dt} = v \\ & \frac{dv}{dt} = u - R v^2 \\ & x(t=0) = 0, ~~ x(t=t_f) = L \\ & v(t=0) = 0, ~~ v(t=t_f) = 0 \\ & -3 \leq u \leq 1 \end{align*}

(1)

where $u$ is the acceleration/braking (your control variable) and $R$ is the drag coefficient (parameter).

Scale Time¶

Let $t = \tau \cdot t_f$ where $\tau \in [0,1]$ . Thus $dt = t_f d\tau$ . The optimal control problem becomes:

\begin{align*} \min_{u,t_f} \quad & t_f \\ \mathrm{s.t.} \quad & \frac{dx}{d\tau} = t_f v \\ & \frac{dv}{d\tau} = t_f (u - R v^2) \\ & x(\tau = 0) = 0, ~~ x(\tau = 1) = L \\ & v(\tau = 0) = 0, ~~ v(\tau = 1) = 0 \\ & -3 \leq u \leq 1 \end{align*}

(2)

Declaring our Model with Pyomo.DAE¶

We can use Pyomo.dae to automatically formulate the collocation equations (i.e., add constraints that numerically integrate the ODE model).

def create_model1():

    # Define the  model
    m = pyo.ConcreteModel()

    # Define the model parameters
    m.R = pyo.Param(initialize=0.001, units=1 / units.m)  # Friction factor
    m.L = pyo.Param(initialize=100.0, units=units.m)  # Final position

    # Define time
    m.tau = dae.ContinuousSet(bounds=(0, 1))  # Dimensionless time set
    m.tf = pyo.Var(initialize=1, units=units.s)  # Final time

    # Define remaining algebraic variables
    m.x = pyo.Var(m.tau, bounds=(0, m.L + 50 * units.m), units=units.m)  # Position
    m.v = pyo.Var(m.tau, bounds=(0, None), units=units.m / units.s)  # Velocity
    m.u = pyo.Var(
        m.tau, bounds=(-3.0, 1.0), initialize=0, units=units.m / units.s / units.s
    )  # Acceleration

    # Define derivative variables.
    #
    # Units MUST be given explicitly: DerivativeVar defaults to dimensionless and
    # does NOT inherit units from the state variable. Because tau is dimensionless,
    # d/dtau has the SAME units as the state itself.
    m.dx = dae.DerivativeVar(m.x, units=units.m)
    m.dv = dae.DerivativeVar(m.v, units=units.m / units.s)

    # Declare the objective (minimize final time)
    m.obj = pyo.Objective(expr=m.tf)

    # Define the constraints
    # position
    def _ode1(m, i):
        if i == 0:
            return pyo.Constraint.Skip
        return m.dx[i] == m.tf * m.v[i]

    m.ode1 = pyo.Constraint(m.tau, rule=_ode1)

    # velocity
    def _ode2(m, i):
        if i == 0:
            return pyo.Constraint.Skip
        return m.dv[i] == m.tf * (m.u[i] - m.R * m.v[i] ** 2)

    m.ode2 = pyo.Constraint(m.tau, rule=_ode2)

    # Define the initial/boundary conditions
    def _init(m):
        yield m.x[0] == 0
        yield m.x[1] == m.L
        yield m.v[0] == 0
        yield m.v[1] == 0

    m.initcon = pyo.ConstraintList(rule=_init)

    # Check dimensional consistency HERE, on the continuous model, BEFORE the
    # discretization equations are added. See ./units_and_pyomo_dae.md and
    # https://github.com/Pyomo/pyomo/issues/1790
    assert_units_consistent(m)

    return m


m = create_model1()

Now let’s inspect the model.

m.pprint()

2 Param Declarations
    L : Size=1, Index=None, Domain=Any, Default=None, Mutable=True, Units=m
        Key  : Value
        None : 100.0
    R : Size=1, Index=None, Domain=Any, Default=None, Mutable=True, Units=1/m
        Key  : Value
        None : 0.001

4 Var Declarations
    tf : Size=1, Index=None, Units=s
        Key  : Lower : Value : Upper : Fixed : Stale : Domain
        None :  None :     1 :  None : False : False :  Reals
    u : Size=2, Index=tau, Units=m/s**2
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          0 :  -3.0 :     0 :   1.0 : False : False :  Reals
          1 :  -3.0 :     0 :   1.0 : False : False :  Reals
    v : Size=2, Index=tau, Units=m/s
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          0 :     0 :  None :  None : False :  True :  Reals
          1 :     0 :  None :  None : False :  True :  Reals
    x : Size=2, Index=tau, Units=m
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          0 :     0 :  None : 150.0 : False :  True :  Reals
          1 :     0 :  None : 150.0 : False :  True :  Reals

1 Objective Declarations
    obj : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : minimize :         tf

3 Constraint Declarations
    initcon : Size=4, Index={1, 2, 3, 4}, Active=True
        Key : Lower : Body : Upper : Active
          1 :   0.0 : x[0] :   0.0 :   True
          2 :     L : x[1] :     L :   True
          3 :   0.0 : v[0] :   0.0 :   True
          4 :   0.0 : v[1] :   0.0 :   True
    ode1 : Size=1, Index=tau, Active=True
        Key : Lower : Body            : Upper : Active
          1 :   0.0 : dx[1] - tf*v[1] :   0.0 :   True
    ode2 : Size=1, Index=tau, Active=True
        Key : Lower : Body                          : Upper : Active
          1 :   0.0 : dv[1] - tf*(u[1] - R*v[1]**2) :   0.0 :   True

1 ContinuousSet Declarations
    tau : Size=1, Index=None, Ordered=Sorted
        Key  : Dimen : Domain : Size : Members
        None :     1 : [0..1] :    2 : {0, 1}

2 DerivativeVar Declarations
    dv : Size=2, Index=tau, Units=m/s
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          0 :  None :  None :  None : False :  True :  Reals
          1 :  None :  None :  None : False :  True :  Reals
    dx : Size=2, Index=tau, Units=m
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          0 :  None :  None :  None : False :  True :  Reals
          1 :  None :  None :  None : False :  True :  Reals

13 Declarations: R L tau tf x v u dx dv obj ode1 ode2 initcon

Aside: checking units, and when to check them¶

Notice where assert_units_consistent(m) sits in create_model1 above: at the very end of the model-building function, before any discretization is applied.

Two details in the code above are worth pausing on.

DerivativeVar does not inherit units. It defaults to dimensionless. A model can look thoroughly annotated and still fail the very first check. Before this notebook declared units= on m.dx and m.dv, assert_units_consistent failed immediately with dimensionless not compatible with meter -- even though every Var and Param had units.
The units of $dx/d\tau$ are metres, not metres per second. We scaled time, so $\tau$ is dimensionless and the derivative carries the same units as the state.

That second point has a pleasant consequence, explored in the cell below.

See Units in Pyomo.dae Models for the full story, including what the check does not look at (Var bounds and .fix() values are never checked).

# The time-scaling trick has a side benefit: because tau is DIMENSIONLESS, the
# finite-difference/collocation coefficients really are dimensionless, so this model
# passes the units check BEFORE and AFTER discretization.
#
# Do not generalize from this! It is a happy accident of scaling time to [0, 1].
# A model whose ContinuousSet is time-in-seconds fails after discretization --
# see notebooks/3/PyomoDAE_TCLab.ipynb for that case.

m_demo = create_model1()  # assert_units_consistent already ran inside
print("Before discretization: units are consistent.")

pyo.TransformationFactory("dae.collocation").apply_to(
    m_demo, nfe=15, scheme="LAGRANGE-RADAU", ncp=3
)

try:
    assert_units_consistent(m_demo)
    print("After discretization:  units are STILL consistent (tau is dimensionless).")
except Exception as e:
    print("After discretization:  FAILED ->", e)

Before discretization: units are consistent.
After discretization:  units are STILL consistent (tau is dimensionless).

Discretize/Transcribe and Solve¶

# discretizer = TransformationFactory('dae.finite_difference')
# discretizer.apply_to(m,nfe=15,scheme='BACKWARD')

# Declare the discretizer
discretizer = pyo.TransformationFactory("dae.collocation")
discretizer.apply_to(m, nfe=15, scheme="LAGRANGE-RADAU", ncp=3)
# discretizer.apply_to(m,nfe=15,scheme='LAGRANGE-LEGENDRE',ncp=3)

# force piecewise constant controls (acceleration) over each finite element
m = discretizer.reduce_collocation_points(m, var=m.u, ncp=1, contset=m.tau)

# Solve
solver = pyo.SolverFactory("ipopt")
results = solver.solve(m, tee=True)
assert pyo.check_optimal_termination(results), (
    f"Solve failed: status={results.solver.status}, "
    f"termination={results.solver.termination_condition}"
)

print("final time = %6.2f seconds" % (pyo.value(m.tf)))

Ipopt 3.14.19: 

******************************************************************************
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.19, running with linear solver MUMPS 5.8.2.

Number of nonzeros in equality constraint Jacobian...:      829
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:      135

Total number of variables............................:      228
                     variables with only lower bounds:       46
                variables with lower and upper bounds:       91
                     variables with only upper bounds:        0
Total number of equality constraints.................:      214
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  1.0000000e+00 1.00e+02 1.30e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1r 1.0000000e+00 1.00e+02 9.99e+02   2.0 0.00e+00    -  0.00e+00 3.26e-07R  8
   2r 1.2449252e+00 9.95e+01 9.93e+02   2.0 1.26e+02    -  5.63e-03 6.22e-03f  1
   3r 2.6155110e+00 9.91e+01 1.22e+02   0.6 1.56e+00    -  7.62e-01 8.77e-01f  1
   4r 5.7407161e+00 9.81e+01 1.16e+02  -0.1 4.77e+01   0.0 2.43e-02 6.55e-02f  1
   5r 6.9131266e+00 9.76e+01 1.19e+02  -0.1 4.77e+01    -  7.69e-02 2.46e-02f  1
   6r 1.8926907e+01 9.25e+01 4.60e+02  -0.1 8.87e+01    -  2.39e-02 1.35e-01f  1
   7r 1.9252497e+01 8.70e+01 3.95e+02  -0.1 1.08e+01    -  2.88e-02 8.03e-01f  1
   8  2.0135727e+01 8.62e+01 6.54e+00  -1.0 2.64e+02    -  1.14e-01 9.09e-03h  5
   9  2.0857378e+01 8.51e+01 1.12e+01  -1.0 3.10e+02    -  1.01e-01 1.20e-02h  6
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  2.1432701e+01 8.41e+01 1.85e+01  -1.0 3.11e+02    -  1.50e-01 1.20e-02h  6

  11  2.2030985e+01 8.29e+01 3.02e+01  -1.0 2.93e+02    -  2.34e-01 1.41e-02h  6
  12  2.2781040e+01 8.14e+01 4.72e+01  -1.0 2.74e+02    -  3.37e-01 1.87e-02h  6
  13  1.4199237e+01 1.24e+02 3.11e+01  -1.0 2.62e+02    -  5.20e-01 4.38e-01H  1
  14  2.1365014e+01 8.49e+01 1.03e+01  -1.0 9.32e+01    -  5.48e-02 2.34e-01h  2

  15  2.2706495e+01 7.92e+01 8.10e+01  -1.0 1.15e+02    -  6.41e-01 6.36e-02h  4
  16  2.3519561e+01 7.41e+01 1.19e+02  -1.0 8.53e+01    -  3.68e-01 6.19e-02h  5
  17  2.3548373e+01 7.38e+01 2.26e+02  -1.0 6.83e+01    -  7.29e-01 3.34e-03h  9
  18  2.0378731e+01 6.54e+01 1.94e+02  -1.0 2.92e+02    -  1.04e-01 1.28e-01h  1
  19  1.9424173e+01 4.16e+01 7.83e+01  -1.0 1.12e+02    -  1.16e-01 3.70e-01h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  20  1.7710108e+01 1.68e+01 3.14e+01  -1.0 9.79e+01    -  3.38e-01 5.75e-01f  1
  21  1.7993344e+01 8.23e+00 2.96e+02  -1.0 2.83e+01    -  1.00e+00 5.14e-01h  1
  22  1.9147366e+01 4.90e-01 1.61e+01  -1.0 9.69e+00    -  1.00e+00 9.90e-01h  1
  23  1.9205734e+01 1.67e-02 1.26e+03  -1.0 4.50e+00    -  1.00e+00 9.90e-01h  1
  24  1.9211008e+01 4.04e-05 1.02e-06  -1.0 1.18e-01    -  1.00e+00 1.00e+00h  1
  25  1.7317872e+01 3.20e+00 1.46e+06  -5.7 1.68e+01    -  7.50e-01 1.00e+00f  1

  26  1.6683635e+01 7.56e-01 1.79e+05  -5.7 1.22e+01    -  8.77e-01 9.88e-01h  1
  27  1.6527179e+01 7.58e-02 3.73e+04  -5.7 7.89e+00    -  7.92e-01 9.92e-01h  1
  28  1.6520375e+01 6.03e-05 5.38e-07  -5.7 8.83e-02    -  1.00e+00 1.00e+00h  1
  29  1.6520269e+01 6.27e-09 8.71e-11  -5.7 9.10e-04    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  30  1.6520192e+01 9.76e-09 1.62e-10  -8.6 1.19e-03    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 30

                                   (scaled)                 (unscaled)
Objective...............:   1.6520191521080950e+01    1.6520191521080950e+01
Dual infeasibility......:   1.6219085743281554e-10    1.6219085743281554e-10
Constraint violation....:   9.7603276572044706e-09    9.7603276572044706e-09
Variable bound violation:   5.4130100402005610e-09    5.4130100402005610e-09
Complementarity.........:   2.5852522114913472e-09    2.5852522114913472e-09
Overall NLP error.......:   9.7603276572044706e-09    9.7603276572044706e-09


Number of objective function evaluations             = 96
Number of objective gradient evaluations             = 26
Number of equality constraint evaluations            = 96
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 32
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 30
Total seconds in IPOPT                               = 0.242

EXIT: Optimal Solution Found.

final time =  16.52 seconds

Plot Results¶

def extract_plot_model1_results(m):

    # Define empty lists
    x = []  # position, units of length
    v = []  # velocity, units of length per time
    u = []  # acceleration, units of length per time squared
    time = []  # time
    tf = pyo.value(m.tf)

    # Loop over time and append the solution values for each variable to their respective lists
    for i in m.tau:
        time.append(i * tf)
        x.append(pyo.value(m.x[i]))
        v.append(pyo.value(m.v[i]))
        u.append(pyo.value(m.u[i]))

    # Make a figure
    plt.figure(figsize=(12, 4))

    # Format subplot 1 (position)
    plt.subplot(131)
    plt.plot(time, x, linewidth=3, label="x")
    plt.title("location (m)", fontsize=16, fontweight="bold")
    plt.xlabel("time (s)", fontsize=16, fontweight="bold")
    plt.tick_params(direction="in", labelsize=15)

    # Format subplot 2 (velocity)
    plt.subplot(132)
    plt.plot(time, v, linewidth=3, label="v")
    plt.xlabel("time (s)", fontsize=16, fontweight="bold")
    plt.title("velocity (m/s)", fontsize=16, fontweight="bold")
    plt.tick_params(direction="in", labelsize=15)

    # Format subplot 3 (acceleration)
    plt.subplot(133)
    plt.plot(time, u, linewidth=3, label="u")
    plt.xlabel("time (s)", fontsize=16, fontweight="bold")
    plt.title("acceleration (m/s$^2$)", fontsize=16, fontweight="bold")
    plt.tick_params(direction="in", labelsize=15)

    plt.show()


extract_plot_model1_results(m)

Discussion Questions

Do the results make sense?

What do you know about the relationships between position, velocity, and acceleration? Do you see this relationship in the plots?

Finer Time Discretization¶

Does the intermediate braking force (see $u \approx -1.2$ m/s/s around $t=12$ s above) go away if we consider a finer time discretization?

m = create_model1()

# Declare the discretizer
discretizer = pyo.TransformationFactory("dae.collocation")
discretizer.apply_to(m, nfe=200, scheme="LAGRANGE-RADAU", ncp=3)

# force piecewise constant controls (acceleration) over each finite element
m = discretizer.reduce_collocation_points(m, var=m.u, ncp=1, contset=m.tau)

# Solve
solver = pyo.SolverFactory("ipopt")
results = solver.solve(m, tee=True)
assert pyo.check_optimal_termination(results), (
    f"Solve failed: status={results.solver.status}, "
    f"termination={results.solver.termination_condition}"
)

print("final time = %6.2f seconds" % (pyo.value(m.tf)))

extract_plot_model1_results(m)

Ipopt 3.14.19: 

******************************************************************************
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.19, running with linear solver MUMPS 5.8.2.

Number of nonzeros in equality constraint Jacobian...:    11004
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:     1800

Total number of variables............................:     3003
                     variables with only lower bounds:      601
                variables with lower and upper bounds:     1201
                     variables with only upper bounds:        0
Total number of equality constraints.................:     2804
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  1.0000000e+00 1.00e+02 2.77e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  1.5275277e+01 9.98e+01 1.02e+02  -1.0 9.82e+03    -  9.91e-05 1.45e-03f  1
   2  3.5437364e+01 8.51e+01 1.14e+05  -1.0 1.38e+02   2.0 1.74e-05 1.48e-01f  1
   3  4.1813834e+01 7.61e+01 9.88e+04  -1.0 1.33e+02   1.5 7.52e-02 1.06e-01f  1
   4  4.1825565e+01 7.61e+01 9.88e+04  -1.0 2.21e+02   1.0 1.39e-01 1.68e-04h  1
   5  9.1089621e+01 6.96e+01 1.38e+05  -1.0 1.17e+02   1.5 2.29e-01 9.90e-01f  1
   6  9.1067800e+01 6.95e+01 1.10e+06  -1.0 2.87e+01   2.8 6.30e-01 7.61e-04h  1
   7  6.4034186e+01 2.10e+01 6.40e+05  -1.0 3.84e+01   2.3 9.56e-01 7.05e-01h  1
   8  6.3919314e+01 2.09e+01 6.20e+05  -1.0 2.32e+01   1.8 9.51e-01 4.95e-03h  1
   9  6.3915285e+01 2.09e+01 5.98e+05  -1.0 1.79e+01   2.3 1.00e+00 7.54e-04h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  5.7448459e+01 9.02e+00 2.38e+05  -1.0 2.62e+01   1.8 6.91e-01 6.45e-01h  1
  11  5.7388561e+01 8.94e+00 2.36e+05  -1.0 1.23e+01   1.3 2.73e-02 9.11e-03h  1

  12  5.7388118e+01 8.93e+00 2.31e+05  -1.0 6.72e+00   1.7 1.00e+00 4.37e-04h  1
  13  5.5324661e+01 6.00e-01 9.44e+03  -1.0 1.07e+01   1.3 6.60e-01 9.91e-01h  1
  14  5.5306135e+01 2.06e-04 1.27e+04  -1.0 1.54e-01   0.8 9.90e-01 1.00e+00f  1
  15  5.5324755e+01 5.39e-05 2.90e-01  -1.0 1.41e-01   0.3 1.00e+00 1.00e+00f  1
  16  5.5308158e+01 2.00e-05 1.72e-02  -3.8 2.51e-02  -0.2 1.00e+00 1.00e+00h  1
  17  5.5251448e+01 2.33e-04 1.90e-02  -3.8 8.30e-02  -0.6 1.00e+00 1.00e+00h  1
  18  5.5158418e+01 7.15e-04 1.51e-01  -3.8 2.40e-01  -1.1 1.00e+00 6.41e-01h  1
  19  5.4800275e+01 9.41e-03 1.85e-02  -3.8 7.30e-01  -1.6 1.00e+00 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  20  5.4630079e+01 1.00e-02 1.13e-01  -3.8 2.01e+00  -2.1 1.00e+00 1.63e-01h  1
  21  5.2283571e+01 4.13e-01 6.29e-02  -3.8 5.59e+00  -2.5 1.00e+00 1.00e+00f  1
  22  5.0371547e+01 6.00e-01 3.26e-01  -3.8 1.48e+01  -3.0 1.00e+00 2.94e-01h  1
  23  4.9652461e+01 6.13e-01 9.04e-01  -3.8 3.52e+01  -3.5 1.00e+00 5.20e-02h  1
  24  4.5660379e+01 1.88e+00 1.26e+00  -3.8 8.65e+01  -4.0 5.64e-01 1.40e-01f  1
  25  4.1365082e+01 3.19e+00 1.24e+00  -3.8 2.34e+02  -4.5 4.20e-04 5.51e-02f  1

  26  4.0197856e+01 3.01e+00 1.17e+00  -3.8 7.93e+01  -4.0 2.28e-02 5.54e-02h  1
  27  3.9964114e+01 2.91e+00 1.13e+00  -3.8 2.47e+01  -3.6 1.00e+00 3.28e-02h  1
  28  3.8325528e+01 2.86e+00 1.04e+00  -3.8 7.88e+01  -4.1 6.16e-02 7.70e-02h  1
  29  3.7141752e+01 2.61e+00 8.65e-01  -3.8 2.39e+01  -3.7 1.00e+00 1.73e-01h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  30  3.6551485e+01 2.52e+00 8.34e-01  -3.8 5.97e+01  -4.1 1.71e-01 3.70e-02h  1
  31  3.5580047e+01 2.07e+00 6.88e-01  -3.8 2.46e+01  -3.7 1.00e+00 1.78e-01h  1
  32  3.4121136e+01 2.16e+00 6.23e-01  -3.8 5.70e+01  -4.2 2.29e-01 9.77e-02h  1
  33  3.3640766e+01 1.99e+00 5.62e-01  -3.8 2.15e+01  -3.8 1.00e+00 9.71e-02h  1
  34  3.2102297e+01 2.00e+00 4.95e-01  -3.8 5.32e+01  -4.2 7.27e-01 1.24e-01h  1
  35  3.0480177e+01 2.20e+00 4.89e-01  -3.8 2.14e+03  -4.7 2.20e-03 1.12e-02h  1
  36  3.0002089e+01 2.13e+00 9.20e+00  -3.8 4.90e+01  -4.3 7.22e-01 4.43e-02h  1

  37  2.9448434e+01 2.14e+00 1.47e+00  -3.8 4.12e+03  -4.8 6.57e-05 2.21e-03h  1
  38  2.8507740e+01 2.11e+00 9.80e+00  -3.8 4.84e+01  -4.3 7.85e-01 9.27e-02h  1
  39  2.8465606e+01 2.08e+00 4.82e+00  -3.8 1.90e+01  -3.9 2.00e-01 1.35e-02h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  40  2.7071467e+01 2.11e+00 3.29e+00  -3.8 4.62e+01  -4.4 1.65e-01 1.59e-01h  1
  41  2.6476677e+01 2.16e+00 3.24e+00  -3.8 9.96e+02  -4.9 8.93e-03 9.46e-03h  1
  42  2.5883551e+01 2.08e+00 1.05e+03  -3.8 4.41e+01  -4.4 1.22e-03 7.61e-02h  1
  43  2.5424859e+01 2.11e+00 1.04e+03  -3.8 1.17e+03  -4.9 5.97e-03 6.83e-03h  1
  44  2.5112008e+01 2.05e+00 9.92e+02  -3.8 4.37e+01  -4.5 1.17e-01 4.31e-02h  1
  45  2.5050205e+01 2.04e+00 9.94e+02  -3.8 2.30e+02  -5.0 3.60e-03 1.92e-03h  1

  46  2.4573469e+01 1.97e+00 9.28e+02  -3.8 4.66e+01  -4.5 5.51e-02 6.57e-02h  1
  47  2.4550117e+01 1.97e+00 9.31e+02  -3.8 2.24e+02  -5.0 2.50e-03 7.44e-04h  1
  48  2.4485029e+01 1.95e+00 9.38e+02  -3.8 7.88e+02  -4.6 1.56e-04 9.53e-03f  1
  49  2.4249827e+01 1.80e+00 1.15e+03  -3.8 2.55e+01  -4.2 8.48e-05 8.43e-02h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  50  2.3949092e+01 1.75e+00 1.10e+03  -3.8 5.21e+01  -4.6 1.85e-01 3.90e-02h  1
  51  2.3651709e+01 1.78e+00 1.10e+03  -3.8 3.17e+03  -5.1 1.55e-03 1.87e-03h  1
  52  2.3084752e+01 1.74e+00 1.02e+03  -3.8 5.27e+01  -4.7 5.11e-02 7.50e-02h  1
  53  2.2814651e+01 1.76e+00 1.02e+03  -3.8 3.50e+03  -5.2 2.98e-03 1.95e-03h  1
  54  2.2683736e+01 1.74e+00 9.97e+02  -3.8 1.93e+02  -4.7 1.32e-03 1.90e-02f  1

  55  2.2327798e+01 1.78e+00 1.01e+03  -3.8 6.88e+03  -5.2 4.61e-05 1.33e-03h  1
  56  2.2020375e+01 1.74e+00 9.58e+02  -3.8 5.47e+01  -4.8 2.95e-03 4.20e-02h  1
  57  2.1985569e+01 1.74e+00 9.57e+02  -3.8 4.45e+02  -5.3 2.93e-02 8.47e-04h  1
  58  2.1433522e+01 1.71e+00 8.79e+02  -3.8 5.62e+01  -4.9 2.23e-01 8.21e-02h  1
  59  2.0739074e+01 1.87e+00 8.51e+02  -3.8 1.95e+02  -5.3 2.52e-02 3.11e-02h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  60  2.0607933e+01 1.83e+00 8.30e+02  -3.8 4.91e+01  -4.9 2.37e-01 2.45e-02h  1
  61  1.9530283e+01 2.33e+00 7.79e+02  -3.8 1.60e+02  -5.4 1.25e-02 6.19e-02h  1
  62  1.9477823e+01 2.30e+00 7.68e+02  -3.8 4.14e+01  -5.0 1.58e-01 1.37e-02h  1
  63  1.8745183e+01 2.54e+00 7.21e+02  -3.8 1.23e+02  -5.4 6.05e-02 6.15e-02h  1
  64  1.8545772e+01 2.37e+00 6.62e+02  -3.8 5.75e+01  -5.0 1.01e-01 8.11e-02h  1

  65  1.8089921e+01 2.42e+00 6.20e+02  -3.8 9.91e+01  -5.5 1.90e-01 6.36e-02h  1
  66  1.7411997e+01 2.95e+00 5.98e+02  -3.8 3.21e+02  -6.0 2.72e-02 3.60e-02h  1
  67  1.7282639e+01 2.99e+00 5.97e+02  -3.8 5.35e+03  -6.4 3.49e-03 1.09e-03h  1
  68  1.7173531e+01 3.00e+00 5.95e+02  -3.8 1.44e+03    -  1.40e-02 3.88e-03h  1
  69  1.6651938e+01 3.60e+00 5.59e+02  -3.8 5.02e+02    -  2.50e-02 6.04e-02h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  70  1.6650651e+01 3.59e+00 5.58e+02  -3.8 4.27e+02    -  1.13e-01 2.55e-03h  1
  71  1.6662814e+01 2.78e+00 4.31e+02  -3.8 7.69e+01    -  2.95e-02 2.26e-01h  1
  72  1.6662429e+01 2.58e+00 4.01e+02  -3.8 4.30e+01    -  2.46e-02 7.09e-02h  1
  73  1.6621587e+01 2.35e+00 3.64e+02  -3.8 1.39e+02    -  3.33e-02 9.07e-02h  1
  74  1.6752920e+01 1.60e+00 2.49e+02  -3.8 4.89e+01    -  1.16e-01 3.17e-01h  1
  75  1.6735208e+01 1.53e+00 2.38e+02  -3.8 1.10e+02    -  9.57e-03 4.51e-02h  1
  76  1.6646743e+01 1.26e+00 1.95e+02  -3.8 7.66e+01    -  3.62e-02 1.77e-01h  1
  77  1.6609125e+01 1.03e+00 1.59e+02  -3.8 4.80e+01    -  2.53e-02 1.86e-01h  1
  78  1.6562492e+01 9.88e-01 1.53e+02  -3.8 2.18e+02    -  5.13e-02 3.89e-02h  1
  79  1.6552168e+01 9.37e-01 1.45e+02  -3.8 3.26e+01    -  4.49e-02 5.12e-02h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  80  1.6490646e+01 8.80e-01 1.36e+02  -3.8 2.15e+02    -  1.65e-02 6.16e-02h  1
  81  1.6499564e+01 8.00e-01 1.24e+02  -3.8 8.76e+00    -  1.00e+00 9.10e-02h  1
  82  1.6558800e+01 1.01e-02 4.90e-03  -3.8 2.76e+00    -  1.00e+00 1.00e+00h  1
  83  1.6558692e+01 2.95e-07 7.10e-09  -3.8 5.58e-02    -  1.00e+00 1.00e+00h  1

  84  1.6473396e+01 5.75e-02 2.41e+03  -5.7 1.11e+01    -  8.38e-01 1.00e+00f  1
  85  1.6471743e+01 6.26e-04 6.34e+01  -5.7 6.23e+00    -  9.74e-01 1.00e+00h  1
  86  1.6471737e+01 2.18e-08 2.02e-11  -5.7 6.29e-02    -  1.00e+00 1.00e+00h  1
  87  1.6470638e+01 1.22e-05 6.67e-01  -8.6 1.83e-01    -  9.96e-01 1.00e+00f  1
  88  1.6470638e+01 2.98e-11 2.55e-14  -8.6 5.63e-04    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 88

                                   (scaled)                 (unscaled)
Objective...............:   1.6470637653941068e+01    1.6470637653941068e+01
Dual infeasibility......:   2.5501024614419560e-14    2.5501024614419560e-14
Constraint violation....:   8.8249407781404443e-12    2.9790692778709322e-11
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   2.5063203258172501e-09    2.5063203258172501e-09
Overall NLP error.......:   2.5063203258172501e-09    2.5063203258172501e-09


Number of objective function evaluations             = 89
Number of objective gradient evaluations             = 89
Number of equality constraint evaluations            = 89
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 89
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 88
Total seconds in IPOPT                               = 1.526

EXIT: Optimal Solution Found.

final time =  16.47 seconds

Yes, the intermediate braking amount gets smaller. The car also reaches the finish line slightly faster.

Another Version¶

Another modeling option is to define an extra differential equation to compute time.

\begin{align*} \min_{u} \quad & t_f \\ \mathrm{s.t.} \quad & \frac{dx}{d\tau} = t_f v \\ & \frac{dv}{d\tau} = t_f (u - R v^2) \\ & \frac{dt}{d\tau} = t_f \\ & x(\tau = 0) = 0, ~~ x(\tau = 1) = L \\ & v(\tau = 0) = 0, ~~ v(\tau = 1) = 0 \\ & t(\tau = 0) = 0 \\ & -3 \leq u \leq 1 \end{align*}

(3)

def create_model2():

    # Define the  model
    m = pyo.ConcreteModel()

    # Define the model parameters
    m.R = pyo.Param(initialize=0.001, units=1 / units.m)  # Friction factor
    m.L = pyo.Param(initialize=100.0, units=units.m)  # Final position

    # Define time
    m.tau = dae.ContinuousSet(bounds=(0, 1))  # Dimensionless time set
    m.time = pyo.Var(m.tau, bounds=(0, None), units=units.s)  # Time
    m.tf = pyo.Var(initialize=1, units=units.s)  # Final time

    # Define remaining algebraic variables
    m.x = pyo.Var(m.tau, bounds=(0, m.L + 50 * units.m), units=units.m)  # Position
    m.v = pyo.Var(m.tau, bounds=(0, None), units=units.m / units.s)  # Velocity
    m.u = pyo.Var(
        m.tau, bounds=(-3.0, 1.0), initialize=0, units=units.m / units.s / units.s
    )  # Acceleration

    # Define derivative variables. tau is dimensionless, so d/dtau carries the
    # same units as the state. DerivativeVar does not infer these -- say them.
    m.dtime = dae.DerivativeVar(m.time, units=units.s)
    m.dx = dae.DerivativeVar(m.x, units=units.m)
    m.dv = dae.DerivativeVar(m.v, units=units.m / units.s)

    # Declare the objective (minimize final time)
    m.obj = pyo.Objective(expr=m.tf)

    # Define the constraints
    # position
    def _ode1(m, i):
        if i == 0:
            return pyo.Constraint.Skip
        return m.dx[i] == m.tf * m.v[i]

    m.ode1 = pyo.Constraint(m.tau, rule=_ode1)

    # velocity
    def _ode2(m, i):
        if i == 0:
            return pyo.Constraint.Skip
        return m.dv[i] == m.tf * (m.u[i] - m.R * m.v[i] ** 2)

    m.ode2 = pyo.Constraint(m.tau, rule=_ode2)

    # time
    def _ode3(m, i):
        if i == 0:
            return pyo.Constraint.Skip
        return m.dtime[i] == m.tf

    m.ode3 = pyo.Constraint(m.tau, rule=_ode3)

    # Define the initial/boundary conditions
    def _init(m):
        yield m.x[0] == 0
        yield m.x[1] == m.L
        yield m.v[0] == 0
        yield m.v[1] == 0
        yield m.time[0] == 0

    m.initcon = pyo.ConstraintList(rule=_init)

    # Units are checked on the continuous model, before discretization.
    assert_units_consistent(m)

    return m


m = create_model2()

# Declare the discretizer
discretizer = pyo.TransformationFactory("dae.collocation")
discretizer.apply_to(m, nfe=15, scheme="LAGRANGE-RADAU", ncp=3)

# force piecewise constant controls (acceleration) over each finite element
m = discretizer.reduce_collocation_points(m, var=m.u, ncp=1, contset=m.tau)

# Solve
solver = pyo.SolverFactory("ipopt")
results = solver.solve(m, tee=True)
assert pyo.check_optimal_termination(results), (
    f"Solve failed: status={results.solver.status}, "
    f"termination={results.solver.termination_condition}"
)

print("final time = %6.2f seconds" % (pyo.value(m.tf)))

Ipopt 3.14.19: 

******************************************************************************
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.19, running with linear solver MUMPS 5.8.2.

Number of nonzeros in equality constraint Jacobian...:     1145
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:      135

Total number of variables............................:      319
                     variables with only lower bounds:       92
                variables with lower and upper bounds:       91
                     variables with only upper bounds:        0
Total number of equality constraints.................:      305
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  1.0000000e+00 1.00e+02 3.48e-01  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  9.9965381e-01 9.99e+01 1.09e+02  -1.0 1.00e+03    -  7.69e-05 1.26e-03f  1
   2  1.0672606e+00 9.98e+01 1.21e+02  -1.0 8.13e+02    -  1.78e-04 1.13e-03f  1
   3  1.6627950e+00 9.94e+01 3.76e+01  -1.0 6.32e+02    -  4.66e-04 3.48e-03f  1
   4  2.1580706e+00 9.89e+01 2.13e+01  -1.0 3.70e+02    -  5.94e-03 4.73e-03f  1
   5  3.3740442e+00 9.75e+01 2.13e+01  -1.0 2.83e+02    -  6.48e-03 1.47e-02f  1

   6  3.4568065e+00 9.75e+01 2.12e+01  -1.0 3.65e+03  -2.0 9.05e-05 2.25e-04f  4
   7  3.5166768e+00 9.75e+01 2.12e+01  -1.0 1.45e+05  -1.6 3.05e-07 2.70e-06f  4
   8  3.6179009e+00 9.74e+01 2.11e+01  -1.0 4.56e+02  -1.1 6.90e-03 7.57e-04h  4
   9  3.7394769e+00 9.74e+01 2.10e+01  -1.0 2.14e+03  -1.6 4.01e-04 2.70e-04h  3
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  3.9133700e+00 9.72e+01 2.10e+01  -1.0 3.70e+02  -1.2 1.66e-02 1.41e-03h  3
  11  4.4185485e+00 9.71e+01 1.98e+01  -1.0 1.58e+03    -  2.59e-03 1.55e-03h  1
  12  4.5952693e+00 9.69e+01 1.98e+01  -1.0 4.16e+02    -  5.05e-03 1.57e-03h  1
  13  4.9840028e+00 9.66e+01 1.97e+01  -1.0 4.67e+02    -  7.60e-03 3.38e-03h  1
  14  5.3631837e+00 9.61e+01 1.98e+01  -1.0 2.80e+02    -  6.53e-03 5.27e-03h  1
  15  5.8588385e+00 9.55e+01 1.96e+01  -1.0 3.63e+02    -  7.97e-03 5.72e-03h  2
  16  6.2740527e+00 9.51e+01 1.96e+01  -1.0 4.21e+02    -  1.56e-02 4.53e-03h  1
  17  6.6359047e+00 9.45e+01 2.20e+01  -1.0 2.85e+02    -  1.32e-02 6.28e-03h  1
  18  7.1860006e+00 9.36e+01 2.75e+01  -1.0 3.17e+02    -  1.47e-02 9.08e-03h  4
  19  7.8876904e+00 9.27e+01 2.76e+01  -1.0 4.06e+02    -  1.29e-02 9.78e-03h  2
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  20  8.3302999e+00 9.20e+01 2.73e+01  -1.0 3.73e+02    -  7.91e-03 7.33e-03h  2
  21  8.7005071e+00 9.15e+01 3.38e+01  -1.0 3.73e+02    -  1.67e-03 5.55e-03h  1
  22  8.8614541e+00 9.12e+01 3.54e+01  -1.0 3.10e+02    -  1.44e-02 3.62e-03h  1
  23  9.3754873e+00 9.01e+01 3.94e+01  -1.0 3.09e+02    -  1.59e-02 1.18e-02h  1
  24  1.0173832e+01 8.84e+01 3.85e+01  -1.0 3.05e+02    -  4.92e-02 1.97e-02h  3

  25  1.1806583e+01 8.59e+01 5.17e+01  -1.0 3.68e+02  -1.7 4.14e-03 2.80e-02h  2
  26  1.3077594e+01 8.20e+01 3.75e+01  -1.0 2.75e+02  -1.3 6.22e-02 4.47e-02h  3
  27  1.3909957e+01 8.03e+01 2.75e+01  -1.0 3.28e+02  -1.7 6.11e-02 2.12e-02h  3
  28  1.4861686e+01 7.65e+01 3.46e+01  -1.0 2.40e+02  -1.3 2.30e-02 4.78e-02f  3

  29  1.5847144e+01 7.35e+01 3.03e+01  -1.0 2.70e+02  -1.8 9.50e-02 3.85e-02h  5
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  30  1.7199412e+01 7.22e+01 2.99e+01  -1.0 5.15e+02  -2.3 2.07e-02 1.84e-02h  1
  31  1.7264467e+01 7.20e+01 4.19e+01  -1.0 2.86e+02    -  1.06e-01 2.56e-03h  1
  32  3.3947489e+01 5.88e+00 1.98e+02  -1.0 2.70e+02    -  7.07e-02 9.22e-01H  1
  33  3.3892282e+01 7.53e-02 7.75e+01  -1.0 1.50e+01  -1.8 5.55e-01 9.90e-01h  1
  34  3.1099781e+01 3.25e+00 4.34e+02  -1.0 2.09e+02    -  1.62e-01 1.83e-01f  1
  35  3.3257932e+01 1.99e+00 9.45e+03  -1.0 2.87e+01    -  5.56e-01 9.92e-01f  1
  36  3.0676688e+01 3.10e+00 1.40e+06  -1.0 3.75e+01    -  2.66e-01 9.12e-01f  1
  37  2.8971974e+01 1.94e+00 5.27e+06  -1.0 3.39e+01    -  2.61e-01 1.00e+00f  1
  38  2.7979354e+01 1.41e+00 9.55e+05  -1.0 2.59e+01    -  7.48e-01 5.84e-01h  1
  39  2.5272864e+01 5.18e+00 1.26e+06  -1.0 1.17e+02    -  1.27e-01 5.83e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  40  2.5693827e+01 3.44e-01 1.27e+06  -1.0 2.11e+01    -  2.36e-01 1.00e+00f  1
  41  2.2982818e+01 4.74e+00 7.78e+05  -1.0 9.26e+01    -  3.86e-01 6.26e-01f  1
  42  2.0165927e+01 6.73e+00 5.59e+05  -1.0 6.11e+01    -  2.81e-01 8.52e-01f  1
  43  1.8548643e+01 5.04e+00 3.42e+05  -1.0 6.52e+01    -  3.89e-01 5.46e-01f  1
  44  1.9691536e+01 1.01e+00 4.17e-02  -1.0 1.57e+01    -  1.00e+00 1.00e+00f  1
  45  1.9758985e+01 5.69e-02 3.33e-03  -1.0 1.62e+01    -  1.00e+00 1.00e+00h  1
  46  1.7116764e+01 5.67e+00 2.62e+06  -2.5 2.17e+01    -  7.30e-01 1.00e+00f  1
  47  1.6751093e+01 3.80e-01 9.60e+04  -2.5 1.33e+01    -  9.63e-01 1.00e+00h  1

  48  1.6637572e+01 4.98e-02 8.50e-05  -2.5 7.31e+00    -  1.00e+00 1.00e+00h  1
  49  1.6637598e+01 7.09e-07 2.95e-08  -2.5 1.37e-01    -  1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  50  1.6521647e+01 2.19e-02 1.17e+04  -5.7 1.83e+00    -  9.59e-01 1.00e+00f  1
  51  1.6520272e+01 3.39e-06 4.17e-08  -5.7 2.75e-02    -  1.00e+00 1.00e+00h  1
  52  1.6520192e+01 9.39e-09 1.64e-10  -8.6 1.01e-03    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 52

                                   (scaled)                 (unscaled)
Objective...............:   1.6520191520761479e+01    1.6520191520761479e+01
Dual infeasibility......:   1.6444820809545745e-10    1.6444820809545745e-10
Constraint violation....:   9.3908738563186489e-09    9.3908738563186489e-09
Variable bound violation:   5.4521387404804500e-09    5.4521387404804500e-09
Complementarity.........:   2.5872224352029854e-09    2.5872224352029854e-09
Overall NLP error.......:   9.3908738563186489e-09    9.3908738563186489e-09


Number of objective function evaluations             = 110
Number of objective gradient evaluations             = 53
Number of equality constraint evaluations            = 110
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 53
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 52
Total seconds in IPOPT                               = 0.445

EXIT: Optimal Solution Found.

final time =  16.52 seconds