Nonlinear model predictive control (NMPC) can be applied to large nonlinear systems with complex objectives while enforcing constraints, but the cost of solving the underlying optimization problems limits real-time application for systems with fast dynamics or many states. Neural networks can reduce this online cost by approximating the control law, but may violate constraints and typically require many precomputed optimizer solutions for training. The learning-based iterative solver for constrained optimization (LISCO) instead learns to predict primal-dual steps toward satisfaction of the Karush-Kuhn-Tucker (KKT) conditions through self-supervised training without precomputed solutions, but its dense, monolithic architecture grows with the prediction horizon. We remove this horizon dependence by exploiting the temporal structure of the problem in the solver architecture, which makes the network size independent of the horizon while leaving the training loss and the solver iteration unchanged. A closed-loop resampling strategy further aligns the training distribution with closed-loop operation at no additional sampling cost. On an economic NMPC benchmark of a CSTR cascade, the proposed solver scales to larger problems than monolithic LISCO and achieves median speedups of up to a factor of four over IPOPT. In closed loop, it achieves comparable control performance to IPOPT with no observed constraint violations.
Figures & tables
Figure 1: Monolithic LISCO update map. From the parameters p and the current iterate z^k , the modified KKT residual ( 8 ) is evaluated and a single dense network maps it to the full primal-dual step Δz^k , which updates the iterate. Iterating this map drives the residual to the tolerance. The decomposed map proposed in this work is shown for comparison in Figure 2 .
Figure 2: LISCO update map with temporal decomposition. The global modified KKT residual ( 8 ) is split into overlapping stage windows. A single shared network is evaluated on all windows in parallel. The overlapping predictions are combined by the fixed averaging operator. The per-stage magnitude and the global step scaling γ are restored to form the primal-dual step.
Figure 3: Time decomposition of the horizon. The stage-partitioned primal-dual iterate z=(z0,…,zN−1) is covered by overlapping windows S0,…,SM−1 of m consecutive stages, shifted by the stride s , so consecutive windows overlap in m−s stages (red). A stage that lies in two windows, such as stage 2 , receives one prediction per window, and these are averaged into a single step by the merge operation ( 25 ). The number of windows is M=1+(N−m)/s .
Symbol
Value [Unit]
Description
nr
study dependent [–]
Number of reactors
Δt
1.0 [s]
Discretization step
VR
15 [L]
Total reactor volume
Vout
1 [L/s]
Volumetric flow rate
k1
2 [L/(mol s)]
Rate constant
k2
2 [L/(mol s)]
Rate constant
Table 1: Parameters of the CSTR cascade model. The first reaction A+B→R is denoted by subscript 1 and the second reaction 2A→S by subscript 2 .
Variable
Unit
Scaling range
Bounds
cA,n
mol/L
[0,0.1]
[0,∞)
cB,n
mol/L
[0,1]
unconstrained
cR,n
mol/L
[0,4]
unconstrained
cS,n
mol/L
[0,1]
unconstrained
Tr,n
∘ C
[20,80]
[20,80]
uA,n
mol/s
[0,1.5]
[0,1.5]
Table 2: Scaling ranges and decision variable bounds. The listed scaling ranges are mapped to [0,1] .
Group
Setting
Value
Network
hidden layers
2×1024 , GELU [ 19 ]
precision
bfloat16
Solver
step scaling γ
0.1
norm clamp ηmax
100
smoothing ϵ
10−1
penalization ρFB
0.8
Table 3: Shared configuration of all studies.
Figure 4: CPU solve time against the horizon length N at nr=5 for the decomposed solver and IPOPT with MA27, with the median (solid) and the maximum (dotted) over 256 instances per configuration. The monolithic solver is not shown, as it does not converge at nr=5 .
monolithic
decomposed
nr
N
npar [M]
ttrain [h]
conv.
npar [M]
ttrain [h]
conv.
2
10
1.84
0.11
1.000
1.47
0.14
1.000
2
20
2.62
0.23
0.945
1.47
0.27
1.000
2
50
4.95
0.54
0.000
1.47
0.77
1.000
2
100
–
–
–
1.47
2.05
1.000
3
10
2.23
0.12
0.996
1.67
0.15
1.000
Table 4: Training cost and convergence of the monolithic and the decomposed solver at equal hidden architecture and training budget for all reactor counts nr and horizon lengths N . Reported are the trainable parameters npar , the training times ttrain and the fraction of the 256 evaluation samples reaching ∥rϕ∥∞≤10−5 within 300 iterations, where 0.000 means that no sample reached the tolerance. Dashes mark configurations without a monolithic run.
tLISCO [ms]
tIPOPT [ms]
speedup [ × ]
nr
N
nz
med
p95
max
med
p95
max
min
med
max
2
10
380
0.89
1.29
1.66
2.65
2.92
3.39
1.64
2.99
4.14
2
20
760
1.69
1.99
6.21
3.47
3.77
4.42
0.52
2.03
2.83
2
50
1900
3.44
3.93
4.36
5.79
6.24
6.99
1.28
1.68
2.23
2
100
3800
6.24
9.91
10.66
10.97
11.63
12.00
1.00
1.72
2.02
3
10
570
1.08
1.46
2.05
3.50
3.81
4.43
1.65
3.23
4.86
Table 5: CPU solve times of the decomposed solver (LISCO) and IPOPT with MA27 over 256 instances per configuration, with the primal-dual dimension nz and the median, 95th percentile and maximum solve times. The speedup is the per-sample ratio tIPOPT/tLISCO on identical instances, reported by its minimum, median and maximum. Both solvers converge on every instance. The monolithic solver is omitted, as it does not converge at most configurations.
Figure 5: Fraction of the 256 samples reaching the residual tolerance ∥rϕ∥∞≤10−5 over the solver iterations for the smoothing parameters ϵ∈{10−1,…,10−6} at nr=2 and N=50 , with an iteration cap of 300.
Figure 6: Absolute relative objective gap ∣q−q∗∣/∣q∗∣ at the final iterate for the smoothing sweep at nr=2 and N=50 over 256 samples, where q∗ denotes the objective of the IPOPT reference solution at tolerance 10−8 . The line marks the median and the boxes span the 10th to 90th percentile.
Figure 7: 90th percentile of the iterations to the residual tolerance ∥rϕ∥∞≤10−5 for different window lengths m and strides s , with nr=5 reactors and N∈{49,50,51} . Each configuration is evaluated on 512 samples with a maximum number of 200 iterations. The grey cell indicates that fewer than 90% of the samples converge within the iteration limit. Dashes indicate combinations that violate the admissibility condition ( 18 ) and were not evaluated. The black outline marks the lowest observed 90th percentile.
metric
IPOPT
LISCO
approx. MPC
closed-loop cost (mean)
-358.4415
-358.8662
-358.4793
solve time median [ms]
18.61
5.20
0.08
solve time p90 [ms]
19.20
6.30
0.10
solve time max [ms]
21.41
11.62
0.43
input viol. (mean)
0.0000
0.0000
0.0000
input viol. (max)
0.0000
0.0000
0.0002
Table 6: Closed-loop comparison of IPOPT, the decomposed solver (LISCO) trained with closed-loop resampling and approximate MPC trained on closed-loop data, over 50 trajectories of 50 steps each, started from uniformly sampled initial states at nr=5 and N=50 . Reported are the mean closed-loop cost, the per-solve times and the mean and maximum constraint violations on the scaled variables for the input and state bounds, where violations below the tolerance 10−5 are not counted. The last two rows list the offline data sampling and training times.
step bin
closed-loop trained
uniform trained
steps 1–1
100.0%
100.0%
steps 2–5
100.0%
84.5%
steps 6–10
100.0%
0.4%
steps 11–25
100.0%
0.0%
steps 26–50
100.0%
0.0%
all 1–50
100.0%
8.8%
Table 7: Fraction of closed-loop solves reaching the residual tolerance ∥rϕ∥∞≤10−5 within 300 iterations, binned by the simulation step, for the decomposed solver trained with closed-loop resampling ( qreset=0.02 ) and with uniform resampling ( qreset=1 ), over the 50 trajectories of 50 steps of the closed-loop study.
Institute of Computer Engineering, TU Wien, Vienna, Austria · AIT Austrian Institute of Technology, Vienna, Austria · Institute of Mechanics and Mechatronics, TU Wien, Vienna, Austria