A Tutorial on Physics Constrained Machine Learning for Hybrid Modeling, Optimization, and Control
Authors: Rahul Golder, Bimol Nath Roy, MM Faruque Hasan
Organizations: Artie McFerrin Department of Chemical Engineering, Texas A&M University, College Station, TX 77843-3122, USA · Texas A&M Energy Institute, Texas A&M University, College Station, TX 77843, USA
Physics-constrained machine learning (PCML) integrates domain knowledge, first-principles, and physical constraints into AI/ML models. In this work, we present a step-by-step tutorial of a modular PCML framework, KKT-HardNet, and illustrate how it could be used for hybrid modeling, process simulation, optimization and model predictive control in a single platform. KKT-HardNet combines fast inference capabilities of neural networks with guaranteed satisfaction of algebraic constraints through differential projection. It can handle heterogeneous data and enforce both linear and nonlinear equality and inequality constraints. The modular design supports custom neural-network backbones and user-defined projection mechanisms for PCML. To enable efficient end-to-end training, KKT-HardNet exploits implicit function theorem to backpropagate through the projection layer. This avoids computationally expensive automatic differentiation. This work also introduces an exact reduced-order Fischer-Burmeister reformulation and a primal-dual interior-point method as an alternative projection solve strategy. We also provide a diverse set of benchmarking problems for PCML in engineering applications. KKT-HardNet achieves orders of magnitude better performance in terms of both constraint satisfaction and inference time. The KKT-HardNet code and the benchmark datasets are available at https://github.com/souls-tamu/kkt-hardnet .
Figures & tables
Figure 1: KKT-HardNet Architecture: The architecture takes parameters as input and through a backbone network predicts the raw output. The projection layer projects the raw output on the constrained surface which creates a feasible output. The model is trained with respective loss function on the constrained output.
Figure 2: Schematic of a neural network.
Figure 3: Graphical interpretation of the projection problem. The unconstrained prediction is pushed towards the feasible domain. While the global minima will represent the minimum distance solution, solving through KKT system guarantees only feasibility. [ 10 ]
Loss function
Type
Application Areas
Ldata=2N1i=1∑N∥y~i−yˉi∥22
Supervised
Hybrid modeling, data-driven optimization, and parameter estimation
Lobj=N1i=1∑Nf(xi,y~i)
Unsupervised
Process optimization, process control
Table 1: Representative loss functions, learning types, and corresponding applications in the KKT-HardNet framework.
Figure 4: Modular organization of KKT-HardNet . Each row corresponds to a stage with different modules.
Figure 5: General KKT-HardNet representation for hybrid modeling, where available process data are combined with partial physical knowledge to learn a constrained input–output mapping.
Figure 7Figure 8
Figure 8: A two-phase Reactor-Condensor Network with recycle
Variable
Description
Nominal value
a
mass transfer area per unit liquid holdup volume (m2/m3)
1000
kA,kC
overall mass transfer coefficients for A and C (mol/(m2⋅s))
20,30
V1T,V2T
total volume of reactor and condenser (m3)
3.0,3.0
ρ
liquid molar density (kmol/m3)
15.0
Table 4: Parameter values and steady-state values for the reactor-condenser system
MSE
Mean eq.
Max eq.
Mean ineq.
Max ineq.
Training Time (s)
Inference Time (ms/batch)
MLP
9.75×10−2
0.29
7.429
2.22×10−5
6.73×10−2
39.78
0.35
PINN
1.31×10−1
0.61
5.29
7.03×10−3
9.80×10−1
41.00
0.31
ENFORCE
7.21×100
7.36×10−10
5.10×10−9
0.00
0.00
16362.92
208.39
KKT-HardNet
3.11×10−3
1.33×10−7
8.59×10−7
6.19×10−8
9.45×10−7
786.21
57.86
Table 5: Performance comparison for the reactor–condenser problem. We define the metrics in appendix A .
Figure 9: General KKT-HardNet formulation for real-time optimization, where the model learns the parametric solution map from the problem parameters θ to a feasible near-optimal solution x~ .
T0=400 K
τ=60 s
kA=5000 /s
kB=106 /s
EA=1×104 cal/mol
EB=1.5×104 cal/mol
R=1.987 cal/(mol K)
ΔH=−5000 cal/mol
ρ=1 kg/L
CP=1000 cal/(kg K)
CA0=1 mol/L
V=100 L
CAs=426.743 K
Qs=40386 cal/s
Table 6: Parameter values and steady-state values for the CSTR case study.
Objective
Mean eq.
Max eq.
Mean ineq.
Max ineq.
Training Time (s)
Inference Time (ms/batch)
IPOPT
0.535
3.43×10−12
9.98×10−11
0.00
0.00
−−
1.83∗
PINN
−0.487
1.96×10−2
4.06×10−2
0.06
0.514
56.87
0.086
DC3
0.560
1.11×10−16
5.71×10−18
0.00
0.00
333.38
2.13
KKT-HardNet
0.535
4.04×10−8
8.57×10−8
0.00
0.00
3725.66
1.65
Table 7: Performance comparison for the CSTR real-time optimization problem in terms of objective value, constraint satisfaction, training time, and inference time for a batch of 100 instances.
Figure 10: Learned RTO policy and its local sensitivity to heat price.
Figure 11: General KKT-HardNet formulation for model predictive control, where the current state, reference, and controller parameters are mapped to the constrained finite-horizon state and control trajectories.
Objective
Mean eq.
Max eq.
Mean ineq.
Max ineq.
Training Time (s)
Inference Time (ms/batch)
IPOPT
1.32×10−1
4.75×10−10
5.63×10−8
4.03×10−14
9.98×10−9
−−
6.30∗
PINN
1.19×10−1
5.31×10−3
8.16×10−2
6.18×10−6
2.72×10−2
137.68
0.72
DC3
2.70×10−1
0.00
0.00
1.39×103
3.45×104
1821.28
31.60
KKT-HardNet
1.60×10−1
2.28×10−5
1.48×10−3
0.00
0.00
10364.44
166.07
Table 8: Performance comparison for the CSTR MPC problem in terms of objective value, constraint satisfaction, training time, and inference time for a batch of 100 instances.
Figure 12: Closed-loop RTO–MPC architecture for the nonlinear CSTR. The RTO layer provides the economically optimal operating target, while the MPC regulates the process toward the updated set point subject to operating and Lyapunov-based constraints.
Figure 13: Closed-loop RTO–MPC response under the time-varying heat-price profile. The RTO layer updates the economically optimal set points (CAsp,CBsp,Tsp,Qsp) , while the KKT-HardNet MPC regulates the nonlinear CSTR toward the updated operating conditions.
Variable
RMSE
Mean settling time (s)
Max. overshoot (%)
Mean steady-state abs. error
CA
6.35×10−4
222.8
0.026
5.27×10−7
CB
6.34×10−4
223.3
0.014
2.85×10−7
T
0.278
204.4
0.004
3.47×10−5
Table 9: Closed-loop tracking performance of the KKT-HardNet MPC controller. Settling time is defined as the time required to enter and remain within 10−4 mol/L of the concentration set points and 0.05 K of the temperature set point until the subsequent RTO update.
Figure 14: Closed-loop Lyapunov analysis of the KKT-HardNet MPC controller. The Lyapunov function increases after each RTO set-point update because the reference changes and subsequently decays toward zero as the process tracks the new target. The Lyapunov margin remains non-negative up to numerical tolerance.
Figure 15: Closed-loop economic performance of the coupled RTO–MPC framework. The top panel compares the instantaneous plant objective Jplant with the steady-state RTO target JRTO , while the remaining panels show the instantaneous and cumulative economic gaps.
Figure 16: A syngas production pathway
MSE
Mean eq.
Max eq.
Mean ineq.
Max ineq.
Training Time (s)
Inference Time (ms/batch)
MLP
6.7×10−5
5.40×10−3
0.11
–
–
92.02
0.2
PINN
1.31×10−4
6.47×10−3
0.05
–
–
90.87
0.2
ENFORCE
3.45×10−5
1.31×10−7
3.81×10−6
–
–
226.70
0.66
KKT-HardNet
4.48×10−5
2.25×10−16
3.55×10−15
–
–
60.80
1.73
Table 10: Performance comparison for the syngas production surrogate-modeling problem in terms of prediction accuracy, elemental-balance constraint satisfaction, training time, and inference time for a batch of 100 instances.
Figure 17: Multi-start projected Adam for syngas operating-condition optimization. More details about the algorithm is given in Appendix E
Figure 18: Optimal syngas operating conditions and corresponding performance metrics as a function of the scalarization parameter α . The left axis shows the scalarized objective, total utility, and CO2 conversion, while the three right-side axes show the optimal feed ratio a , operating pressure Pop , and operating temperature Top .
Appendix figures & tables3 assets
Supplementary material from the paper’s appendix.
Appendix
Figure 19: LMPC control-effort analysis during closed-loop RTO–LMPC operation. The first panel compares the applied heat input Q with the RTO steady-state target Qsp . The second panel shows the LMPC correction u0=Q−Qsp , which is largest immediately following an RTO set-point update and subsequently decays toward zero. The third panel shows the change between successive implemented heat inputs, with ΔQmax=5×103 cal/s shown as a reference. The final panel shows the cumulative absolute LMPC control effort, ∫∣u0∣dt .
Profile
CA RMSE
T RMSE (K)
Max. ∣T−Tsp∣ (K)
Max. ∣u0∣ (cal/s)
Integrated positive economic gap
Max. plant violation
Baseline
6.35×10−4
0.278
1.972
13.81
0.702
0.00
Sinusoidal
1.10×10−3
0.761
6.345
45.66
1.451
0.00
Random
3.87×10−3
1.845
20.216
139.31
4.238
0.00
Spike
3.41×10−3
1.580
22.053
163.46
2.149
0.00
Ramp
7.90×10−4
0.299
1.880
13.16
0.000
0.00
Appendix
Table 11: Closed-loop performance under alternative heat-price profiles.
Figure 20: Closed-loop RTO–LMPC performance under challenging heat-price profiles. The first panel shows the imposed baseline, sinusoidal, random, spike, and ramp heat-price trajectories. The subsequent panels show the corresponding CA tracking error, temperature tracking error, LMPC correction u0 , and positive economic gap. More abrupt price changes generate larger transient tracking errors and corrective actions, while all cases remain bounded and satisfy the physical plant limits.