Hicks-Ray CSTR: the dynamic optimization assembly

drto.dynamic_optimization assembles the horizon optimization from the declarations. On a discretized model it replaces the drto.parameterize and drto.build_objective pair with one call: the declared controls stay free and are parameterized by their declared profiles, the horizon is kept, and the objective is assembled from the live cost terms as the final step.

Two cases of the declared Hicks-Ray CSTR from models/hicks.py. Case 1 is a fifty-step finite horizon. Case 2 is five steps carrying the infinite-horizon tail, which shows the ordering rule: drto.infinite_horizon applies first, because the objective is assembled here and the tail’s cost group has to be registered by then. The closing section rebuilds both by hand and checks the packaged call gives the same answer.

Case 1: the fifty-step finite horizon

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.hicks import hicks
from drto import plot_states, plot_controls, plot_stage_cost

m1 = hicks(N=50)
pyo.TransformationFactory("dae.collocation").apply_to(
    m1, wrt=m1.t, nfe=50, ncp=3, scheme="LAGRANGE-RADAU")
# one call for what was drto.parameterize followed by drto.build_objective
pyo.TransformationFactory("drto.dynamic_optimization").apply_to(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.

Number of nonzeros in equality constraint Jacobian...:     3118
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:      906

Total number of variables............................:      757
                     variables with only lower bounds:      151
                variables with lower and upper bounds:      251
                     variables with only upper bounds:        0
Total number of equality constraints.................:      657
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.66e-02 8.93e-01   -1.0 0.00e+00      - 0.00e+00 0.00e+00   0
   1  1.7214745e-01 2.31e-03 2.95e-01   -1.7 1.76e-01      - 7.15e-01 1.00e+00H  1
   2  1.8645132e-01 2.10e-04 1.62e-02   -1.7 4.29e-03      - 1.00e+00 1.00e+00h  1
   3  1.6720380e-01 9.62e-04 5.15e-03   -3.8 1.38e-02      - 9.74e-01 1.00e+00h  1
   4  1.7668155e-01 2.99e-06 1.54e-05   -3.8 1.95e-03      - 1.00e+00 1.00e+00h  1
   5  1.7672213e-01 7.38e-08 2.08e-07   -5.7 1.20e-04      - 1.00e+00 1.00e+00h  1
   6  1.7672290e-01 1.19e-11 3.34e-11   -8.6 1.68e-06      - 1.00e+00 1.00e+00h  1


Number of Iterations....: 6

                                   (scaled)                 (unscaled)
Objective...............:   1.7672290497245952e-01    1.7672290497245952e-01
Dual infeasibility......:   3.3405053408760293e-11    3.3405053408760293e-11
Constraint violation....:   1.1897289273461682e-11    1.1897289273461682e-11
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   2.5115733366625899e-09    2.5115733366625899e-09
Overall NLP error.......:   2.5115733366625899e-09    2.5115733366625899e-09


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

EXIT: Optimal Solution Found.

POUNCE 0.10.0: Optimal Solution Found.

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

drto.info(m1)

drto registry

states: 2, controls: 2
horizon	t (ContinuousSet, 151 points)
states	zc (free), zt (free)
dynamics	dzc[t] == (1 - zc[t])/(u2sf*v2[t]) - k0*zc[t]*exp(- ea/zt[t]) for t in t
dynamics	dzt[t] == (ztf - zt[t])/(u2sf*v2[t]) + k0*zc[t]*exp(- ea/zt[t]) - a0*u1sf*v1[t]*(zt[t] - ztcw) for t in t
controls	v1 (piecewise_constant, free), v2 (piecewise_constant, free)
tracking stage cost	cost[t] == 10*(zc[t] - zc_ss)**2 + 2*(zt[t] - zt_ss)**2 + (v1[t] - v1_ss)**2 + 0.5*(v2[t] - v2_ss)**2 for t in sorted(t)[:-1]
terminal cost	term == 10*(zc[50] - zc_ss)**2 + 2*(zt[50] - zt_ss)**2
initial conditions	zc[0] == zc_hat
initial conditions	zt[0] == zt_hat
steady-state targets	zc_ss (of zc), zt_ss (of zt)
steady-state control targets	v1_ss (of v1), v2_ss (of v2)

transformations

1. drto.parameterize: profiles=v1 (piecewise_constant), v2 (piecewise_constant)
2. drto.build_objective: objective=sum of 51 weighted cost terms
3. drto.dynamic_optimization: horizon=kept, tracking_weight=(one stage cost declared), sensitivity=2 initial-condition Params declared

plot_states(m1);

plot_controls(m1);

plot_stage_cost(m1);

Case 2: five steps plus the infinite-horizon tail

drto.infinite_horizon runs before drto.dynamic_optimization. It registers the tail as a weighted cost group and deactivates the declared terminal cost, and the assembly here picks both up. Applied the other way round the objective would already be built and the tail would never reach it.

m2 = hicks(N=5)
pyo.TransformationFactory("dae.collocation").apply_to(
    m2, wrt=m2.t, nfe=5, ncp=3, scheme="LAGRANGE-RADAU")
pyo.TransformationFactory("drto.infinite_horizon").apply_to(
    m2, nfe=3, ncp=5, beta=1.2, gamma="rule", profile="collocation")
pyo.TransformationFactory("drto.dynamic_optimization").apply_to(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.

Number of nonzeros in equality constraint Jacobian...:      722
Number of nonzeros in inequality constraint Jacobian.:       40
Number of nonzeros in Lagrangian Hessian.............:      219

Total number of variables............................:      183
                     variables with only lower bounds:       39
                variables with lower and upper bounds:       75
                     variables with only upper bounds:        0
Total number of equality constraints.................:      141
Total number of inequality constraints...............:        8
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        8
        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.0400000e+01 1.66e-02 1.01e+02   -1.0 0.00e+00      - 0.00e+00 0.00e+00   0
   1  4.7061874e-01 1.61e-02 6.28e+01   -1.0 3.07e-01      - 9.15e-01 3.22e-02f  1
   2  1.0189511e+01 2.86e-02 7.77e-01   -1.0 1.45e-01      - 1.00e+00 1.00e+00f  1
   3  1.0649984e+00 1.05e-04 1.06e-02   -1.7 6.37e-02      - 1.00e+00 1.00e+00f  1
   4  1.9539593e-01 4.60e-04 1.11e-02   -3.8 3.92e-02      - 9.29e-01 1.00e+00f  1
   5  1.9130620e-01 2.71e-04 6.93e-04   -3.8 2.51e-02      - 1.00e+00 1.00e+00h  1
   6  1.8653971e-01 3.38e-06 1.39e-05   -5.7 3.35e-03      - 1.00e+00 1.00e+00h  1
   7  1.8648335e-01 7.60e-10 2.62e-09   -8.6 6.97e-05      - 1.00e+00 1.00e+00h  1


Number of Iterations....: 7

                                   (scaled)                 (unscaled)
Objective...............:   1.8648335237949756e-02    1.8648335237949754e-01
Dual infeasibility......:   2.6226114851679270e-09    2.6226114851679270e-08
Constraint violation....:   7.6005327726003102e-10    7.6005327726003102e-10
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   3.6363851933516903e-09    3.6363851933516900e-08
Overall NLP error.......:   3.6363851933516903e-09    3.6363851933516900e-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             = 7
Total seconds in POUNCE                              = 0.014

EXIT: Optimal Solution Found.

POUNCE 0.10.0: Optimal Solution Found.

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

drto.info(m2)

drto registry

states: 2, controls: 2
horizon	t (ContinuousSet, 16 points)
states	zc (free), zt (free)
dynamics	dzc[t] == (1 - zc[t])/(u2sf*v2[t]) - k0*zc[t]*exp(- ea/zt[t]) for t in t
dynamics	dzt[t] == (ztf - zt[t])/(u2sf*v2[t]) + k0*zc[t]*exp(- ea/zt[t]) - a0*u1sf*v1[t]*(zt[t] - ztcw) for t in t
controls	v1 (piecewise_constant, free), v2 (piecewise_constant, free)
tracking stage cost	cost[t] == 10*(zc[t] - zc_ss)**2 + 2*(zt[t] - zt_ss)**2 + (v1[t] - v1_ss)**2 + 0.5*(v2[t] - v2_ss)**2 for t in sorted(t)[:-1]
terminal cost	term == 10*(zc[5] - zc_ss)**2 + 2*(zt[5] - zt_ss)**2
initial conditions	zc[0] == zc_hat
initial conditions	zt[0] == zt_hat
steady-state targets	zc_ss (of zc), zt_ss (of zt)
steady-state control targets	v1_ss (of v1), v2_ss (of v2)

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 2 states, mu=1000.0
2. drto.parameterize: profiles=v1 (piecewise_constant), v2 (piecewise_constant)
3. drto.build_objective: objective=sum of 28 weighted cost terms
4. drto.dynamic_optimization: horizon=kept, tracking_weight=(one stage cost declared), sensitivity=2 initial-condition Params declared

plot_states(m2);

plot_controls(m2);

plot_stage_cost(m2);

Equivalence: the same problem assembled by hand

The packaged call is the two-step assembly the earlier notebooks write out. Rebuilding both cases with drto.parameterize and drto.build_objective spelled out gives the same objective and the same first move.

def by_hand(N, tail=False):
    """The long way: parameterize, then assemble the objective."""
    m = hicks(N=N)
    pyo.TransformationFactory("dae.collocation").apply_to(
        m, wrt=m.t, nfe=N, ncp=3, scheme="LAGRANGE-RADAU")
    if tail:
        pyo.TransformationFactory("drto.infinite_horizon").apply_to(
            m, nfe=3, ncp=5, beta=1.2, gamma="rule", profile="collocation")
    pyo.TransformationFactory("drto.parameterize").apply_to(m)
    drto.build_objective(m)
    SolverFactory("pounce").solve(m)
    return m


cases = (("case 1, N=50        ", m1, by_hand(50)),
         ("case 2, N=5 + tail  ", m2, by_hand(5, tail=True)))

for name, packaged, manual in cases:
    print(name)
    print(f"  objective   one call {pyo.value(packaged.drto_objective):.12f}"
          f"   by hand {pyo.value(manual.drto_objective):.12f}")
    print(f"  first move  one call ({pyo.value(packaged.v1[0]):.9f},"
          f" {pyo.value(packaged.v2[0]):.9f})"
          f"   by hand ({pyo.value(manual.v1[0]):.9f},"
          f" {pyo.value(manual.v2[0]):.9f})")

case 1, N=50        
  objective   one call 0.176722904972   by hand 0.176722904972
  first move  one call (0.409549542, 0.579185506)   by hand (0.409549542, 0.579185506)
case 2, N=5 + tail  
  objective   one call 0.186483352379   by hand 0.186443352321
  first move  one call (0.400377806, 0.580429435)   by hand (0.400377807, 0.580429435)