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.
While physics-informed neural networks (PINNs) have shown strong potential for process modeling, physical equations are only enforced as soft constraints during training, and thus, they do not guarantee constraint satisfaction at inference. We propose a framework, called piecewise-linear Karush--Kuhn--Tucker hard-constrained PINNs (PL-KKT-hPINNs), that strictly enforces nonlinear equality constraints through piecewise-linear projection. This extends the KKT-hPINN framewor, which exactly enforces linear equalities through the Karush--Kuhn--Tucker (KKT) conditions associated with orthogonally projecting neural network outputs onto the constraint feasible region. The method is demonstrated on a continuous stirred-tank reactor (CSTR) case study for both one and two inputs. Results show that PL-KKT-hPINN preserves predictive accuracy comparable to that of a standard neural network while achieving substantially lower constraint violations. In addition, the proposed model shows improved robustness in low-data regimes, yielding lower RMSE than the unconstrained neural network for limited training sample sizes. These results demonstrate that PL-KKT-hPINN provides a computationally efficient and physically consistent framework for surrogate modeling of nonlinear chemical engineering systems.
Fateme Mohammad Mohammadi, Hector Budman, Joshua L. Pulsipher
Department of Chemical Engineering University of Waterloo, 200 University Ave W, Waterloo, ON N2L 3G1, Canada
Hybrid models that combine physics-based and data-driven components have shown strong potential for achieving accuracy and interpretability in control applications. While recent methods have made progress in incorporating physical consistency, challenges remain in scalability, robustness to noise, and control of model complexity. This paper proposes a Physics-Encoded Modular Hybrid Layer (PE-MHL) framework, in which a baseline physics-based model is incrementally refined through the addition of new sub-models, where each new component adds complexity while preserving what previous components have already learned. We establish a theoretical guarantee for this construction: with a least-squares initialization of each new sub-model, the training error is monotonically non-increasing in the number of sub-models and provably converges. Empirical evaluations on a nonlinear NARX benchmark and the Quanser Aero 2 platform demonstrate that PE-MHL matches or outperforms its fully connected monolithic relaxation (over five times more parameters) in accuracy and generalization, while also providing more stable training dynamics and better preservation of underlying data structures.
Ismail Hassaballa, Mircea Lazar
Eindhoven University of Technology, Eindhoven, Netherlands
Generative models have emerged as scalable surrogates for physical simulation, yet they offer no guarantee that their outputs respect the conservation laws, boundary conditions, and nonlinear invariants that govern the underlying physics. Constrained sampling closes this gap, enforcing such constraints exactly at inference time without retraining, but at a computational cost: projection, correction and trajectory-optimization steps are repeated during sampling, with these steps becoming expensive for nonlinear constraints. Standard ML frameworks exacerbate this: their dense tensor algebra and limited sparse solver composability obscure the structure that physical constraints naturally induce, making efficient batched nonlinear optimization difficult to realize in practice. We address this bottleneck by exploiting the structure that sample-wise batching and local PDE couplings induce in the projection subproblems -- namely, block-sparse Jacobian and KKT systems -- exposing this structure using ExaModels.jl and solving the resulting sparse nonlinear programs with MadNLP.jl and GPU sparse factorization. Applied to Physics-Constrained Flow Matching (PCFM), on PDE benchmarks with linear, nonlinear, one-dimensional, and two-dimensional constraints, this approach accelerates nonlinear constraint projection while maintaining constraint satisfaction. These results show that sparse GPU nonlinear optimization is a practical foundation for constrained generative sampling in scientific machine learning.