First-order system: the infinite horizon vs a long finite horizon

The model is the canonical declared first-order system from models/first_order.py: see the module for its construction. Two cases of the same model. Case 1 runs five sampling steps with the infinite-horizon terminal segment (drto.infinite_horizon) carrying the tail, which deactivates the declared terminal cost since the tail owns the cost-to-go. Case 2 runs fifty steps of plain finite horizon with the terminal cost live.

Case 1: five samples plus the infinite-horizon tail

import pyomo.environ as pyo
import pyomo_pounce  # registers pounce with the native factory
from pyomo.contrib.solver.common.factory import SolverFactory
import drto
from models.first_order import first_order
from drto import plot_states, plot_controls, plot_stage_cost

m1 = first_order(N=5)
pyo.TransformationFactory("dae.collocation").apply_to(m1, wrt=m1.t, nfe=5, ncp=3, scheme="LAGRANGE-RADAU")
pyo.TransformationFactory("drto.infinite_horizon").apply_to(
    m1,
    nfe=3,                    # finite elements on the terminal segment
    ncp=5,                    # Gauss-Legendre points per element
    beta=1.2,                 # tail overestimation factor, > 1
    gamma="rule",             # the mesh rule, tanh(gamma*h) = tau_11; pass a number to override
    profile="collocation",    # segment control profile (pyomo-cvp)
)
pyo.TransformationFactory("drto.parameterize").apply_to(m1)  # the finite-horizon profiles
drto.build_objective(m1)
SolverFactory("pounce").solve(m1, tee=True)

********************************************************************************

                    ####    ###   /   # /#   #/  ####  #####
                    #   #  #   # /#   #/ ##  /  #      #
                    ####   #   #/ #   /  # #/#  #      ####
                    #      #   /  #  /#  # /##  #      #
                    #       ##/    #/#   #/  #   ####  #####

********************************************************************************
This program contains POUNCE, a pure-Rust interior-point optimization solver
for nonlinear, conic, and global problems (its NLP core is ported from Ipopt).
Released under the Eclipse Public License (EPL) — drop-in compatible with Ipopt.
         For more information visit https://github.com/jkitchin/pounce
********************************************************************************

This is POUNCE version 0.10.0, running with linear solver FERAL.

Reading C:\Users\Devin\AppData\Local\Temp\tmp0kipog41\unknown1nuk1eh1.nl...
Parsed 90 vars, 73 cons, jac_nnz=325, h_nnz=41 in 0.02s
Problem class: NLP. Selected solver: NLP filter line-search interior-point (pounce-nlp) [solver_selection=auto].

Number of nonzeros in equality constraint Jacobian...:      305
Number of nonzeros in inequality constraint Jacobian.:       20
Number of nonzeros in Lagrangian Hessian.............:       41

Total number of variables............................:       90
                     variables with only lower bounds:        2
                variables with lower and upper bounds:       20
                     variables with only upper bounds:        0
Total number of equality constraints.................:       69
Total number of inequality constraints...............:        4
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        4
        inequality constraints with only upper bounds:        0

iter      objective   inf_pr   inf_du lg(mu)    ||d|| lg(rg) alpha_du alpha_pr  ls
   0  4.4595491e+01 1.00e-01 1.02e+02   -1.0 0.00e+00      - 0.00e+00 0.00e+00   0
   1  2.3009452e+01 9.74e-02 6.89e+01   -1.0 1.89e+00      - 7.79e-01 3.23e-02f  1
   2  5.7603845e+00 1.31e-01 7.92e-01   -1.0 1.14e-01    0.0 1.00e+00 1.00e+00f  1
   3  5.4870343e-01 2.70e-03 3.30e-02   -1.7 1.23e-01      - 1.00e+00 1.00e+00h  1
   4  1.6528570e-01 2.90e-03 2.83e-08   -2.5 4.32e-02      - 1.00e+00 1.00e+00h  1
   5  1.1763246e-01 2.75e-04 1.50e-09   -3.8 1.80e-02      - 1.00e+00 1.00e+00h  1
pounce: wrote C:\Users\Devin\AppData\Local\Temp\tmp0kipog41\unknown1nuk1eh1.sol
   6  1.1524063e-01 1.61e-06 1.85e-11   -5.7 2.27e-03      - 1.00e+00 1.00e+00h  1
   7  1.1520672e-01 7.19e-11 4.26e-14   -8.6 1.53e-05      - 1.00e+00 1.00e+00h  1


Number of Iterations....: 7

                                   (scaled)                 (unscaled)
Objective...............:   1.1520672332537066e-02    1.1520672332537066e-01
Dual infeasibility......:   4.2632564145606011e-14    4.2632564145606011e-13
Constraint violation....:   7.1908492181582684e-11    7.1908492181582684e-11
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   2.5468980039820931e-09    2.5468980039820930e-08
Overall NLP error.......:   2.5468980039820931e-09    2.5468980039820930e-08


Number of objective function evaluations             = 8
Number of objective gradient evaluations             = 8
Number of equality constraint evaluations            = 8
Number of inequality constraint evaluations          = 8
Number of equality constraint Jacobian evaluations   = 8
Number of inequality constraint Jacobian evaluations = 8
Number of Lagrangian Hessian evaluations             = 8
Total seconds in POUNCE                              = 0.010

EXIT: Optimal Solution Found.

POUNCE 0.10.0: Optimal Solution Found.
Status: Solve_Succeeded

<pyomo.contrib.solver.common.results.Results at 0x224a1b46030>

drto.info(m1)

drto registry

states: 1, controls: 1
horizon	t (ContinuousSet, 16 points)
states	z (free)
dynamics	dzdt[t] == (- z[t] + u[t])/tau_p for t in t
controls	u (piecewise_constant, free)
tracking stage cost	cost[t] == 10*(z[t] - z_ss)**2 + (u[t] - u_ss)**2 for t in sorted(t)[:-1]
terminal cost	term == 10*(z[5] - z_ss)**2
initial conditions	z[0] == z_hat
steady-state targets	z_ss (of z)
steady-state control targets	u_ss (of u)

transformations

1. drto.infinite_horizon: segment=3 elements x 5 Legendre points, beta=1.2, gamma=0.01563827, profile=collocation, horizon=kept, infinite tail appended, terminal_cost=terminal deactivated (the tail owns it), terminal=soft pin z(tau=1)=z_s on 1 state, mu=1000.0
2. drto.parameterize: profiles=u (piecewise_constant)
3. drto.build_objective: objective=sum of 24 weighted cost terms

plot_states(m1);

plot_controls(m1);

plot_stage_cost(m1);

Case 2: fifty samples, finite horizon with the terminal cost

m2 = first_order(N=50)
pyo.TransformationFactory("dae.collocation").apply_to(m2, wrt=m2.t, nfe=50, ncp=3, scheme="LAGRANGE-RADAU")
pyo.TransformationFactory("drto.parameterize").apply_to(m2)
drto.build_objective(m2)
SolverFactory("pounce").solve(m2, tee=True)

********************************************************************************

                    ####    ###   /   # /#   #/  ####  #####
                    #   #  #   # /#   #/ ##  /  #      #
                    ####   #   #/ #   /  # #/#  #      ####
                    #      #   /  #  /#  # /##  #      #
                    #       ##/    #/#   #/  #   ####  #####

********************************************************************************
This program contains POUNCE, a pure-Rust interior-point optimization solver
for nonlinear, conic, and global problems (its NLP core is ported from Ipopt).
Released under the Eclipse Public License (EPL) — drop-in compatible with Ipopt.
         For more information visit https://github.com/jkitchin/pounce
********************************************************************************

This is POUNCE version 0.10.0, running with linear solver FERAL.

Reading C:\Users\Devin\AppData\Local\Temp\tmp504a07gn\unknownrdjhzgly.nl...
Parsed 401 vars, 351 cons, jac_nnz=1348, h_nnz=100 in 0.02s
Problem class: NLP. Selected solver: NLP filter line-search interior-point (pounce-nlp) [solver_selection=auto].

Number of nonzeros in equality constraint Jacobian...:     1348
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:      100

Total number of variables............................:      401
                     variables with only lower bounds:        0
                variables with lower and upper bounds:       50
                     variables with only upper bounds:        0
Total number of equality constraints.................:      351
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 1.00e-01 5.00e-01   -1.0 0.00e+00      - 0.00e+00 0.00e+00   0
   1 -1.1178489e+01 2.82e-01 1.67e+00   -1.7 2.33e-01      - 8.42e-01 1.00e+00f  1
   2 -2.0213792e+00 4.60e-02 1.27e-02   -1.7 1.87e-01      - 9.65e-01 1.00e+00h  1
   3  1.1300008e-01 5.49e-05 1.34e-15   -2.5 4.59e-02      - 1.00e+00 1.00e+00h  1
   4  1.1522596e-01 4.10e-07 1.08e-15   -3.8 4.99e-04      - 1.00e+00 1.00e+00h  1
   5  1.1522669e-01 8.13e-10 1.10e-15   -5.7 2.25e-05      - 1.00e+00 1.00e+00h  1
pounce: wrote C:\Users\Devin\AppData\Local\Temp\tmp504a07gn\unknownrdjhzgly.sol
   6  1.1522670e-01 1.03e-13 5.37e-16   -8.6 2.55e-07      - 1.00e+00 1.00e+00h  1


Number of Iterations....: 6

                                   (scaled)                 (unscaled)
Objective...............:   1.1522669525418940e-01    1.1522669525418940e-01
Dual infeasibility......:   5.3739791787211870e-16    5.3739791787211870e-16
Constraint violation....:   1.0277802914293588e-13    1.0277802914293588e-13
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   2.5070789087396480e-09    2.5070789087396480e-09
Overall NLP error.......:   2.5070789087396480e-09    2.5070789087396480e-09


Number of objective function evaluations             = 7
Number of objective gradient evaluations             = 7
Number of equality constraint evaluations            = 7
Number of inequality constraint evaluations          = 7
Number of equality constraint Jacobian evaluations   = 7
Number of inequality constraint Jacobian evaluations = 7
Number of Lagrangian Hessian evaluations             = 7
Total seconds in POUNCE                              = 0.008

EXIT: Optimal Solution Found.

POUNCE 0.10.0: Optimal Solution Found.
Status: Solve_Succeeded

<pyomo.contrib.solver.common.results.Results at 0x224a484f390>

drto.info(m2)

drto registry

states: 1, controls: 1
horizon	t (ContinuousSet, 151 points)
states	z (free)
dynamics	dzdt[t] == (- z[t] + u[t])/tau_p for t in t
controls	u (piecewise_constant, free)
tracking stage cost	cost[t] == 10*(z[t] - z_ss)**2 + (u[t] - u_ss)**2 for t in sorted(t)[:-1]
terminal cost	term == 10*(z[50] - z_ss)**2
initial conditions	z[0] == z_hat
steady-state targets	z_ss (of z)
steady-state control targets	u_ss (of u)

transformations

1. drto.parameterize: profiles=u (piecewise_constant)
2. drto.build_objective: objective=sum of 51 weighted cost terms

plot_states(m2);

plot_controls(m2);

plot_stage_cost(m2);