A Tutorial on Physics Constrained Machine Learning for Hybrid Modeling, Optimization, and Control

arXiv:2610.32914v1 [cs.LG] 23 Oct 2026

A Tutorial on Physics Constrained Machine Learning for Hybrid Modeling, Optimization, and Control

Rahul Golder Artie McFerrin Department of Chemical Engineering Texas A&M University, College Station, TX 77843-3122, USA rahulgolder8420@tamu.edu Bimol Nath Roy Artie McFerrin Department of Chemical Engineering Texas A&M University, College Station, TX 77843-3122, USA bimolnathroy@tamu.edu M. M. Faruque Hasan ∗ 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 hasan@tamu.edu

Abstract 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. § https://github.com/souls-tamu/kkt-hardnet

3 pip install kkt-hardnet

Keywords: Scientific Machine Learning, Hybrid Modeling, Physics Informed Neural Networks, Optimization, Process Control ∗ Corresponding author: hasan@tamu.edu

Preprint. October 26, 2026

1

Contents 1 2

7 A B C D E

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . KKT-HardNet: A Physics Constrained Machine Learning Framework . . . . . . . . . . . . . . . 2.1 Neural Network Backbone . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Projection Layer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Loss Function and Backpropagation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 A Modular Python Package (pip install kkt-hardnet) . . . . . . . . . . . . . . . . . PCML for Hybrid Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Hybrid Modeling of Chemical Processes . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Hybrid Modeling of a Two-Phase Reactor-Condensor System . . . . . . . . . . . . . . . . PCML for Real-Time Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 A CSTR Optimization Case Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . PCML-based Model Predictive Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 A CSTR Control Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . PCML for Grey-Box Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1 Simulation-based Optimization of a Syngas Production Process . . . . . . . . . . . . . . . Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Benchmarking: metrics and other methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Methods Included to Solve the Projection Problem . . . . . . . . . . . . . . . . . . . . . . . . . Reduction and Heuristics on KKT Conditions via Fishcer-Burmeister Reformulation . . . . . . . Additional Closed-Loop Analysis of the RTO–LMPC Framework . . . . . . . . . . . . . . . . . Multi-Start Projected Adam Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

Introduction

3

4 5 6

2 4 4 6 7 8 11 11 12 16 16 19 20 27 27 32 36 36 37 38 40

Mathematical models provide representations of real-world problems that help describe systems behavior and enable the analysis and identification of suitable solutions. However, constructing accurate mathematical models can be difficult for complex systems, particularly when the underlying phenomena and mechanisms are only partially known or computationally expensive to represent. In such cases, data-driven learning offers an attractive alternative, while existing domain knowledge can still provide valuable structure and physical consistency. Purely data-driven machine learning such as multilayer perceptrons (MLPs) and deep neural networks (DNNs) can provide fast approximations, but its performance is often limited by scarce data, poor extrapolation, and the lack of guarantees that predictions remain physically meaningful. Incorporating domain knowledge, such as governing equations, conservation laws or other prior information, can reduce the effective search space and improve data efficiency, robustness, and generalization [1, 2]. Scientific machine learning (SciML) has increasingly explored the integration of physical constraints into neural-networks (NNs). NNs are effective function approximators and have been widely used for surrogate modeling, physics-informed learning [3], inverse design, operator learning, and optimization. However, standard neural-network training is fundamentally unconstrained, which makes the direct enforcement of physical and operational constraints challenging. Penalty-based approaches address this by augmenting the loss function with weighted constraint violations, but fixed penalty weights generally do not guarantee exact feasibility or satisfaction of the stationarity conditions of the corresponding constrained problem. This limitation has motivated the development of physics-constrained machine learning (PCML) methods. Different approaches have been proposed in literature to incorporate physics in machine learning models to some extent. These can be broadly categorized into two classes: (a) soft constrained neural networks and (b) 2

hard constrained neural networks. Soft constrained neural networks, such as physics-informed neural networks (PINNs) [2], incorporate physical constraints by penalizing physics violations through the loss function during training. Due to the simplicity and ability to incorporate differential constraints, this method and its variants have been successfully used to model fluid flow [4, 5, 6], biological systems [7, 8], heat transfer [4], and uncertainty quantification [9], among others, where first-principles are available. This involves balancing between the data-fidelity and the first-principle violation through weights: 𝑁

min 𝛉

1 ∑︁ L (𝐹NN (𝛉, 𝒙)) + 𝜆 𝑅 𝑅 + 𝜆 𝑝 𝐿 physics 𝑁 𝑖=1

(1)

where, 𝐿 physics is the physics loss to penalize the violation of the known first-principles or physical constraints, and 𝜆 𝑝 is the corresponding penalty weight. Soft constrained methods minimize both the objective function (data loss) and the constraint violation (physics loss) simultaneously as an unconstrained function with pre-defined multipliers for each loss term. Thus it diverges from the concept of minimizing the Lagrangian where the multipliers are also variables of the minimization problem. To this point, soft constrained neural networks never guarantee to satisfy the physical constraints. This often leads to unreliable prediction or inference. Unlike soft constrained approaches, hard constrained neural networks aim to enforce constraints strictly by design. As modeling physical systems sometime require exact satisfaction of algebraic constraints that are highly nonlinear in nature, hard constrained methods have gained great attention in recent years in the scientific community. PCML combines the strengths of computationally efficient prediction capability of data driven neural networks with guaranteed satisfaction of physical laws and high accuracy through projection, thereby providing both accurate and efficient hybrid modeling capabilities [10]. These methods aim to overcome the limitations of MLPs and PINNs. For example, systems having safety critical constraints are usually inequality constraints that are not well suited for PINNs. The methodological approaches for hard constrained neural networks can be classified into two general categories: (i) predict-and-complete and (ii) projection-based methods [11]. The predict-and-complete approaches first predict a subset of variables and completes the rest from the equality constraints [12, 13]. An advance to this approach involves completion through explicit or implicit set of equations and then correction through internal gradient descent procedure to enforce inequalities [14]. Projection-based methods exploits optimization theory to restore feasibility. The strategy involves first predicting unconstrained variables and then projecting them to a constrained manifold. Advances in this method include strict satisfaction of linear equalities [15], linear inequalities [16, 17], convex constraints [18], nonlinear equalities and inequalities [10, 19, 20, 21, 22]. In this tutorial we focus on the projection-based method to enforce general nonlinear constraints. There are several key advantages of projection based hard-constrained neural networks. First projectionbased hard-constrained methods separate feasibility restoration from loss minimization. The projection layer enforces the constraints, while the neural-network backbone is trained primarily to minimize the loss, thereby reducing the burden on gradient-based learning. Second, the framework has a nested optimization structure: the inner problem projects the network output onto the feasible set, while the outer problem trains the network to minimize the fitting objective. This structure suggests potential connections to bilevel optimization that warrant further investigation. Third, the approach can extend to parametric optimization, where the active constraint set varies with the problem parameters. Through training across parameter realizations, the network can implicitly learn these active-set transitions, potentially providing useful approximations to the critical-region structure studied in multiparametric programming [23]. With this realization, one can attempt to use hard-constrained neural networks to learn parametric solutions of optimization problems [20], where training encourages the neural network to generate predictions close to both feasible and optimal solution manifold. These predictions can therefore serve as high-quality warm starts for conventional optimization 3

solvers, potentially reducing the number of iterations required for convergence. Hard-constrained learning also provides a promising approach for repeatedly solving parametric optimization problems at low inference cost, with applications in areas such as process control, optimal power flow, and trajectory optimization where rapid online decisions are required. In this article, we present a tutorial of PCML with efficient and modular implementation of KKT-HardNet [10], a neural-network-based PCML framework that enforces constraint feasibility through projection. We provide a step by step tutorial to allow users on how to use the KKT-HardNet code and the framework for various process system engineering applications. Among the novel contributions, we extend the original framework with multiple projection strategies, including the reduced-order Fischer–Burmeister formulation and a primal–dual interior-point method. To enable efficient end-to-end training, gradients through the converged projection are computed using implicit differentiation instead of autodifferentiation through the full sequence of solver iterations. We further provide an open-source, modular implementation with support for sparse computations to exploit special structure in large-scale problems. Finally, for deployment scenarios dominated by single-instance predictions, we support ahead-of-time (AOT) compilation to reduce runtime overhead and enable low-latency inference. The tutorial is organised as follows: Section 2 provides a foundational description of KKT-HardNet and its components along with implementation details. With the basic idea of the PCML framework, we then demonstrate how it can be used for hybrid modeling of complex system 3, real-time and online optimization 4, model predictive control 5, and PCML embedded greybox optimization 6. Across these applications, KKTHardNet is evaluated in terms of predictive accuracy, constraint satisfaction, and computational performance, with comparisons against existing state-of-the-art PCML approaches.

2

KKT-HardNet: A Physics Constrained Machine Learning Framework

In this section we describe a projection-based PCML framework, KKT-HardNet along with brief overview of different components. The framework is provided in Figure 1. Given an unconstrained neural-network prediction, the projection layer solves an euclidian projection problem to obtain a nearby feasible point. This architecture can retain prediction accuracy while achieving constraint violations close to machine precision. We begin with exploring neural networks and their types followed by a brief introduction to projection layer and the training algorithm. Though machine learning covers a broad variety of tools, we limit our discussion in the space of neural networks. Here, we briefly describe the different components of KKT-HardNet: a NN backbone, projection layer and the loss function. Note that, KKT-HardNet is a general framework where different architectures of neural network and projection algorithms can be integrated. 2.1

Neural Network Backbone

Neural networks (NNs) are a model structure to capture complex and nonlinear systems, originally inspired from mimicking how brains process information [24]. Despite being originated from neuroscience, NNs are now used widely due to their impressive capability of universally approximating any arbitrary function [25] (universal approximation theorem, see [26] for details). NNs are perceived as networks consist of artificial neurons (known as perceptrons) arranged in layers [27]. The simplest structured NNs are MLPs shown in Figure 2a. Mathematically, neurons are made with a composition of linear and nonlinear functions, making it suitable for approximating general nonlinear functions. The layers in NNs are arranged sequentially and the output of each layer goes to the next layer as input. The first and last layers are input and output layers and the intermediate layers are hidden layers. The information given to the input layer are the input features (𝒙 ∈ R𝑛 ) and the information coming out of the output layer are the predicted outputs ( 𝒚ˆ ∈ R𝑚 ), where “hat" represents unconstrained output. Let ℎ 𝑘𝑗 be the output of 𝑗 th neuron in 𝑘 th hidden layer (shown in Figure 2b), which is computed as 4

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.

Input Layer

Hidden Layer 1

Hidden Layer 2

Output Layer

𝑦1

⋮

𝑥1

Hidden Layer K

ℎ1𝑘−1

𝑘 𝑤1𝑗

𝑏𝑗𝑘

𝑥2

𝑦2

ℎ2𝑘−1

𝑘 𝑤2𝑗

⋮

∑

⋮

⋮

⋮ ⋮

⋮

𝑥𝑛

⋮

⋮

ℎ𝑝𝑘−1

𝑦𝑚

𝑘 𝑤𝑝𝑗

(a) Multilayer perceptron

Activation function 𝑘 th layer 𝑗 th neuron 𝑤: weights 𝑏: bias

(b) Neuron

Figure 2: Schematic of a neural network.

5

𝜎

ℎ𝑗𝑘

ℎ 𝑘𝑗 = 𝜎 𝑗𝑘

𝑆∑︁ 𝑘−1

! 𝑤 𝑖𝑘 ℎ𝑖𝑘−1 + 𝑏 𝑘𝑗

(2)

,

𝑖=1

where 𝑆 𝑘−1 is the number of neurons in the (𝑘 − 1) th layer (𝑆 𝑘−1 = 𝑝 in Figure 2b), ℎ𝑖𝑘−1 are the inputs coming from (𝑘 − 1) 𝑡 ℎ layer, 𝑤 𝑖𝑘𝑗 are the weights for the 𝑖 , 𝑗 neuron connection, 𝑏 𝑘𝑗 is the bias associated with the particular neuron, and 𝜎 is the activation function. For a MLP with 𝐾 layers, the functional form of the NN, output is represented as follows:  𝒚 = 𝑾 𝐾+1 𝒉 𝐾 + 𝒃 𝐾+1 ,    𝑘  𝒚ˆ = 𝐹NN (𝛉, 𝒙) = 𝒉 = 𝝈 𝑘 (𝑾 𝑘 𝒉 𝑘−1 + 𝒃 𝑘 ), ∀𝑘 = 1, 2, . . . , 𝐾,    𝒉0 = 𝒙, 

(3)

where 𝛉 represents the NN parameters, i.e., the weights and biases over which the NN is trained with a loss function and 𝒉 𝑘 , 𝝈 𝑘 and 𝒃 𝑘 are the hidden layer outputs, the activation functions, and the biases of the neurons in 𝑘 th layer respectively. The matrix 𝑾 𝑘 are the weights for the neuron connection between 𝑘 th and (𝑘 − 1) th layers. Typically, 𝐾 + 1 refers to the output layer where generally the use of activation function is avoided (except for special cases). Commonly used activation functions are sigmoid, rectified linear unit (ReLU), hyperbolic tangent (tanh), etc. (see [28] for details). The choice of activation function depends on the modeling purpose and the nonlinearity in the data, among others. 2.2

Projection Layer

Before defining the projection layer, we recall that our goal is to train a neural network 𝒚 = 𝐹NN (𝛉, 𝒙) while satisfying the general nonlinear constraints of the form 𝒉(𝒙, 𝒚) = 0 and 𝒈(𝒙, 𝒚) ≤ 0, where 𝒙 ∈ R𝑛 is the neural network input, 𝒚 ∈ R𝑚 is the neural network output, 𝒉 : R𝑛+𝑚 ↦→ R 𝑝 and 𝒈 : R𝑛+𝑚 ↦→ R𝑞 are the set of equality and inequality constraint functions with 𝑝 and 𝑞 constraints respectively. We make the following assumptions: 1. the number of equality constraints is less than or equal to the number of neural network output, i.e., 𝑝 ≤ 𝑚 to ensure sufficient degrees of freedom for learning, and 2. the constraint manifold Sfeasible = {𝑦|𝒉( 𝒙¯ , 𝒚) = 0, 𝒈( 𝒙¯ , 𝒚) ≤ 0} is nonempty for a given 𝒙¯ and holds linear independence constraint qualification (LICQ). ෝ 𝒚

𝒉 𝒙, 𝒚 = 𝟎 𝒈 𝒙, 𝒚 ≤ 𝟎 Feasible Region ෥ 𝒚 𝒉 𝒙, 𝒚 = 𝟎 𝒈 𝒙, 𝒚 ≥ 𝟎

𝒈 𝒙, 𝒚 ≤ 𝟎

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] 6

Let, 𝒚ˆ = 𝐹NN (𝛉, 𝒙¯ ) be the raw output from the neural network for a given 𝒙¯ , the projection layer is defined as the following feasibility restoration step (commonly used in nonlinear optimization [29]): 1 ∥ 𝒚 − 𝒚ˆ ∥ 2 𝒚 2 s.t. ℎ𝑖 ( 𝒙¯ , 𝒚) = 0, 𝑖 = 1, 2, . . . , 𝑝 𝑔 𝑗 ( 𝒙¯ , 𝒚) ≤ 0, 𝑗 = 1, 2, . . . , 𝑞,

𝒚˜ = P( 𝒚ˆ , 𝒙) ¯ = arg min

(4)

where, P( 𝒚ˆ , 𝒙) ¯ is the implicit projection step. Mathematically, the solution to this optimization problem is the minimum distance feasible solution from the unconstrained and generally infeasible point 𝒚ˆ (see Figure 3). 𝑦˜ is the ultimate output that is guaranteed to meet the physical constraints. While analytical solution to the projection problem in Eq. 4 exists for only linear equality constraints, solving the problem for general nonlinear constraints is not trivial and require iterative steps. One way is considering the Karush-Kuhn-Tucker (KKT) conditions of the problem to find a feasible solution (in some cases the optimal solution) to guarantee constraint feasibility. See Appendix B for brief details on solving the projection problem. While the efficiency of the forward projection problem is subject to the solution strategy used, efficiency of training relies on the backpropagation. Implicit function theorem can be utilized to make the training efficient with implicit backward gradient computation. With implicit differentiation, the projection algorithm may contain non-differentiable intermediate operations, provided that the converged solution admits a locally differentiable implicit representation. This substantially broadens the class of optimization algorithms that can be incorporated into the projection layer [30]. 2.3

Loss Function and Backpropagation

While the projection layer ensures that the final prediction, 𝑦˜ , strictly satisfies the equality and inequality constraints, we now present an updated loss function for the NN, based on 𝑦˜ , so that the NN also learns the physics. To that end, while various loss functions exist (see [31] for a broader discussion), training KKT-HardNet boils down to solving the structural risk minimization (SRM) problem of the following form: 𝑁

min 𝜽

1 ∑︁ L 𝜃 𝒚˜ 𝑖 + 𝜆 𝑅 𝑅(𝜽). 𝑁 𝑖=1

(5)

where 𝒚˜ 𝑖 = 𝐹NN (𝜽, 𝒙 𝑖 ) is the model prediction of the 𝑖 𝑡 ℎ sample, 𝑁 is the number of training samples, L denotes the task-specific loss function, 𝑅(𝜽) is a regularization term, and 𝜆 𝑅 is the corresponding regularization weight. The choice of L depends on the learning task. For supervised surrogate modeling and parameter-estimation problems, the loss typically measures the discrepancy between the model prediction and the available target data. For optimization and control applications, the physical or economic objective of the underlying optimization problem can instead be used directly as the training loss. Training, therefore, consists of determining the network parameters 𝜽 that minimize the selected task-specific objective while accounting for regularization. Representative loss functions and their corresponding applications considered in this work are summarized in Table 1. One important point to note that, depending on the choice of loss function in the PCML framework, one could deploy PCML for a wide range of process systems engineering capabilities (some are listed in Table 1), including hybrid model development, learning solutions of parametric optimization problems, estimating parameters, etc, by steering the 𝑦ˆ towards the desired outcome. Finally projection corrects 𝑦ˆ to 𝑦˜ . Training NNs was popularized by the success of the error backpropagation algorithm [32] which is a firstorder algorithm to train deep neural networks. It computes the gradient of the total loss with respect to the NN parameters through chain rule and utilizes the gradient to update the parameters. While backpropagation is efficient in gradient computation, second order methods such as Gauss-Newton [33], Levenberg-Marquardt [34, 35] algorithms have also shown promise to train NNs. 7

Table 1: Representative loss functions, learning types, and corresponding applications in the KKT-HardNet framework. Loss function

Type

Application Areas

Supervised

Hybrid modeling, data-driven optimization, and parameter estimation

Unsupervised

Process optimization, process control

𝑁

1 ∑︁ 2 Ldata = 𝒚˜ 𝑖 − 𝒚¯ 𝑖 2 2𝑁 𝑖=1 𝑁

Lobj =

2.4

1 ∑︁ 𝑓 𝒙 𝑖 , 𝒚˜ 𝑖 𝑁 𝑖=1

A Modular Python Package (pip install kkt-hardnet) STAGE 1: INPUT DATA • Importing data • Input and output • Splitting data

ANN

TARGET • Simulation Model • Parameter estimation • Optimization proxy

SETTINGS

PROJECTION

• ANN specifications • Activation function • Loss function

• Variables • Constraint sets • Solution strategy

• • •

Training setup Projection setup Configurations

STAGE 2: MODEL BUILDER

Assemble components

Build model

Ready to train model

STAGE 3: MODEL TRAINER 𝜑(𝑥)

Initialize

Forward pass

Compute loss

Backpropagation

Update weights

Checkpoint

Trained Model

STAGE 4: RESULTS SAVE RESULTS

FUTURE USE

• Training history • Metrics & plots • Model artifacts

• Load trained model • Generate new prediction • Analyze results

Figure 4: Modular organization of KKT-HardNet. Each row corresponds to a stage with different modules. KKT-HardNet is available as a Python package (pip install kkt-hardnet) and extends the original hard-coded proof-of-concept framework [10] into a modular implementation designed for practical use and extension. The package uses fixed-shape JAX computations so that the neural-network and projection operations can be compiled and executed together on a selected device [36]. Figure 4, shows the four stages of the workflow: input, model building, model training, and results. The input stage specifies the data, constraints, neural-network architecture, and training/projection configurations; the model builder assembles these components; the training stage optimizes the neural-network backbone through the projected outputs while keeping the projection layer non-trainable; and the results stage saves the trained model, artifacts,

8

and training history for future use. To support this workflow, KKT-HardNet provides a unified interface for defining neural-network backbones, differentiable constraints, projection methods, loss functions, and training procedures through six core modules: problems, nn, losses, opt, train, and export. Their responsibilities are summarized below: The major modules of KKT-HardNet package are: • problems: This module includes robust ways of defining physics and operational constraints starting from symbolically representing the variables. Equations can be defined individually and as blocks with matrix form. We provide options for defining specific structured constraints such as linear or quadratic constraints. While simulation/surrogate modeling does not involve any optimization objective, for optimization proxy this module includes an option to define objective. Furthermore, variable bounds can be included separately as a special set of constraints. • nn: This module utilizes jax.nn and includes common functions for NN libraries to define the NN backbone of the model. User has control over the architecture including selection of the activation function and the size of the NN backbone. Some commonly used activation functions are sigmoid, tanh, relu, etc. • losses: This module deals with defining task specific loss function. For example, simulation/surrogate modeling, and parameter estimation involves data fitting for which mean squared error loss is defined with a consistency loss. For optimization, this creates a loss function that incorporates the objective function of the optimization problem. • opt: This module includes different strategies for solving the projection problem. In our implementation, we have three methods: (i) explicit projection when only linear equality constraints are present, (ii) solving KKT system with Fishcer-Burmeister reformulation which gives a KKT-based feasible solution, and (iii) primal-dual interior point algorithm 2 [37]. While the explicit projection is a single step computation after unconstrained prediction, methods (ii) and (iii) involve iterative procedure to find a solution. Backward gradient computation is expensive for large scale problems through automatic differentiation. In out implementation, we leverage implicit function theorem to compute the implicit backward gradient for the projection layer (see [30, 18, 20] for details on implicit backward gradient). We also introduced a reduced method to solve KKT system with Fishcer-Burmeister reformulation which we describe in Appendix C. • train: The training stage takes the assembled pipeline and trains the NN backbone to minimize specified loss function. The ‘train’ module includes the training and validation data along with the user defined training configurations to train and evaluate the model. This module also checks for the best found loss and saves that model for downstream use. • export: This module provides the functionality required to serialize, compile, and load trained KKT-HardNet predictors for deployment. The implementation supports two artifact types: (i) a portable predictor represented using the JAX export format, and (ii) an environment-specific native predictor obtained by compiling the portable artifact for the current hardware and software environment. The portable loader deserializes the exported JAX computation and wraps it using a common predictor interface, allowing the model to be evaluated without reconstructing the original training code. A command line interface is additionally provided to compile a portable predictor into a native artifact directly on the target machine. The native artifact is target specific and requires a compatible JAX and jaxlib runtime. Therefore, it should not be interpreted as a completely standalone operating-system executable. Codebox 1 illustrates the corresponding problem-definition workflow. GPU support is available through pip install "kkt-hardnet[cuda12]". 2 We use the available open source code of the implementation to reproduce in high level python to make it suitable for batch operations and to

include in training loop with implicit backward gradient.

9

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54

% ====================================================== % Step 1: Import Necessary Libraries % ====================================================== import os os.environ["JAX_ENABLE_X64"] = "True" from pathlib import Path import numpy as np import jax jax.config.update("jax_enable_x64", True) import jax.numpy as jnp # Imports from KKT-HardNet from kkthn import KKTHardNet, Loss, MLP, Projection % ====================================================== % Step 2: Define Configurations % ====================================================== SEED = 42 TRAIN = {# Define your Configuration here} DATA_PATH = "your-data.csv" % ====================================================================== % Step 3: Model Components: Backbone, Projection Layer and Loss Function % ====================================================================== # Backbone nn = MLP(hidden_sizes=(64, 64), activation="relu") # Projection Layer projection = Projection(name="my_model", method="fischer-burmeister-reduced") x = projection.add_parameter(["x1", "x2", "x3"]) y = projection.add_variable(["y1", "y2", "y3"]) projection.constraints.add( y.y1 + y.y2 + y.y3 == 1, y.y1 >= 0, y.y1 <= 10, # Define other constraints) # Loss Function loss = Loss.mse() % ====================================================================== % Step 4: Model Definition % ====================================================================== model = KKTHardNet("surrogate", nn, loss, projection, train=TRAIN) model.dataset(parameters=DATA_PATH, variables=DATA_PATH) % ====================================================================== % Step 3: Model Training % ====================================================================== result = model.model() model.summary() model.plot_history()

Code 1: Code Snippet for using KKT-HardNet. In the subsequent sections, we explore different applications of KKT-HardNet across different engineering fields. We include application codes to provide a self-contained tutorial and minimize reliance on external resources. Across all the case studies we present, a significant portion of the code remains same. These include library imports, backbone setup, training configurations, and solver configurations. See appendix A for the details on the performance metrics and the methods we consider to benchmark KKT-HardNet against. 10

3

PCML for Hybrid Modeling

Hybrid modeling combines data-driven learning with available first-principles knowledge to develop predictive models of complex systems. This is particularly useful for engineering processes where a complete mechanistic model may be difficult or computationally expensive to solve, while process data are available from experiments, simulations, or sensors. Although a neural network can learn the input–output relationship directly from these data, its predictions are not inherently guaranteed to satisfy known physical relations. KKT-HardNet provides a systematic way to incorporate such partial knowledge as hard constraints while retaining the flexibility of data-driven learning. In general, let 𝒙 ∈ R𝑛𝑥 denote the observed process inputs and 𝒚¯ ∈ R𝑛𝑦 denote the corresponding measured or simulated outputs. A neural-network backbone learns the mapping 𝒚˜ = 𝐹NN (𝜽, 𝒙) from the available data, while known physical relations are represented by ℎ(𝒙, 𝒚) = 0 and 𝑔(𝒙, 𝒚) ≤ 0. KKT-HardNet projects the backbone prediction onto this constraint set, such that the final prediction 𝒚 satisfies the incorporated physics. Importantly, the complete governing model is not required; only the available physical knowledge needs to be included in the projection layer. This allows the remaining relationships to be learned from data and makes the framework suitable for partially known systems. The general formulation of hybrid modeling with KKT-HardNet is illustrated in 5.

𝒙

System

𝒚

Data ഥ 𝒚

𝐿𝑑𝑎𝑡𝑎

Inputs 𝒙

..

yƶ

Physics 𝒉 𝒙, 𝒚 = 𝟎 𝒈(𝒙, 𝒚) <= 𝟎

𝐿𝑐𝑜𝑛𝑠𝑖𝑠𝑡𝑒𝑛𝑐𝑦

Outputs 𝒚ƿ Satisfies 𝒉 𝒙, 𝒚 = 𝟎 𝒈(𝒙, 𝒚) <= 𝟎

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.

3.1

Hybrid Modeling of Chemical Processes

Here, we show the numerical and computational improvement of the framework compared to the base study reported by Iftakher et al. [10] (available at github.com/souls-tamu/kkt-hardnet/tree/old-workspace). We compare two case studies: (i) Chemical process involving extractive distillation and (ii) Pooling problem involving nonlinear equality and inequality. For details of the systems we refer to [10]. Code for these case studies are available at the GitHub repository. 11

3.1.1

A case study on chemical process simulation involving nonlinear thermodynamics

Key information: 3 input variables, 9 output variables, 6 equality constraints, zero inequality constraints, includes nonlinearity. Table 2: Performance comparison on the extractive distillation column case study

Condenser 𝐹𝐷 𝐷 𝑥𝑅32 𝐷 𝑥𝑅125 𝐷 𝑥𝐼𝐿

IL solvent 𝐿 𝑅𝑅

𝐹𝐼𝐿

KKT-HardNet (this tutorial)

Extractive Distillation Column R-410A

𝐹𝐹 𝐹 𝑥𝑅32 𝐹 𝑥𝑅125

Reboiler

𝐹𝐵 𝐵 𝑥𝑅32 𝐵 𝑥𝑅125 𝐵 𝑥𝐼𝐿

Figure 6: Extractive distillation of R-410A using an Ionic liquid [EMIM][SCN]. 3.1.2

Metric

Base

FB-N

FB-R

PDIP

T. MSE V. MSE T. Viol. V. Viol. Time (s)

1.05 × 10−4 1.08 × 10−4 2.0 × 10−8 2.0 × 10−8 12381.94

1.10 × 10−7 3.16 × 10−8 1.12 × 10−7 5.53 × 10−8 122.92

5.67 × 10−8 3.07 × 10−8 6.60 × 10−7 7.05 × 10−7 98.76

3.27 × 10−8 1.72 × 10−8 7.87 × 10−13 9.71 × 10−13 108.36

Abbreviations FB-N: fischer-burmeister-newton FB-R: fischer-burmeister-reduced PDIP: primal-dual-interior-point T. MSE: Train MSE; V. MSE: Validation MSE T. Viol.: Train Violation; V. Viol.: Validation Violation Violation is measured as the mean absolute violation. Base studies from [10] were run on the same hardware for fair comparison.

A case study on pooling problem involving nonlinear equality and inequality constraints

Key information: 4 input variables, 5 output variables, 4 equality constraints, 2 inequality constraints, includes nonlinearity. Feed 𝐴 3% Sulfur Feed 𝐵 1% Sulfur

𝑃𝑥 Pool

Mixer

Product 𝑋 ≤2.5% Sulfur

𝑚% Sulfur

KKT-HardNet (this tutorial)

𝑃𝑦 𝐶𝑥

Feed C 2% Sulfur

Table 3: Performance comparison on the pooling problem case study

Metric Product 𝑌 ≤1.5% Sulfur

T. MSE V. MSE Figure 7: Process flowsheet of the pool- T. Viol. ing problem. The product can be formed V. Viol. Time (s) from three feeds. 3.2

Splitter

𝐶𝑦

Mixer

Base

FB-N

FB-R

PDIP

49.65 69.11 1.53 × 10−5 1.53 × 10−5 8224.61

54.99 45.43 1.51 × 10−7 1.44 × 10−7 103.72

55.05 43.20 2.55 × 10−7 3.25 × 10−7 76.44

54.98 43.49 4.99 × 10−11 6.53 × 10−11 68.07

Hybrid Modeling of a Two-Phase Reactor-Condensor System

To illustrate the hybrid modeling capability of a complex process simulation, we consider a two-phase reactor–condenser system with liquid recycle [38], (illustrated in Figure 8). The reactor pressure dynamics and interphase mass-transfer dynamics are assumed to be fast. Reactants 𝐴 and 𝐵 are fed to the reactor at molar flow rates 𝐹𝐴 and 𝐹𝐵 and temperatures 𝑇𝐴 and 𝑇𝐵 , respectively, with 𝐴 entering in the vapor phase and 𝐵 in the liquid phase. Inside reactor, reactant 𝐴 transfers to the liquid phase at a rate 𝑁 𝐴1 , where the exothermic reaction 𝐴 + 𝐵 → 2𝐶 occurs.   −𝐸 𝑎 𝑅 𝐴 = 𝑘 10 exp 𝑀1𝑙 𝜌𝑥 𝐴1 𝑥 𝐵1 , (6) 𝑅𝑇1 12

𝐹2𝑣 , 𝑦𝐴2

𝐹1𝑣 𝑀2𝑣 , 𝑦𝐴2 𝑇2 , 𝑃2

𝑀1𝑣 , 𝑦𝐴1 𝑇1 , 𝑃1

𝐹𝐴 , 𝑇𝐴

𝑁𝐴2

𝑁𝐶1

𝑁𝐴1

𝐹2𝑙

𝑁𝐶2

𝑀2𝑙 𝑥𝐴2

𝑄2

𝐹𝐵 , 𝑇𝐵 𝑀1𝑙

𝑥𝐴1 𝑥𝐵1

𝑄1

𝐴 + 𝐵 → 2𝐶

𝐹1𝑙

Figure 8: A two-phase Reactor-Condensor Network with recycle Table 4: Parameter values and steady-state values for the reactor-condenser system Variable

Description

Nominal value

𝑎 𝑘 𝐴, 𝑘 𝐶 𝑉1𝑇 , 𝑉2𝑇 𝜌

mass transfer area per unit liquid holdup volume (𝑚 2 /𝑚 3 ) overall mass transfer coefficients for A and C (𝑚𝑜𝑙/(𝑚 2 · 𝑠)) total volume of reactor and condenser (𝑚 3 ) liquid molar density (𝑘𝑚𝑜𝑙/𝑚 3 )

1000 20, 30 3.0, 3.0 15.0

The training dataset was generated by solving a DAE model at steady state for a range of operating conditions. The problem property values are given in table 4. Six process inputs, 𝐹𝐴 and 𝐹𝐵 , 𝑇𝐴 and 𝑇𝐵 , 𝐹2𝑙 , 𝐹1𝑙 , were varied using Latin hypercube sampling to provide a stratified coverage of the prescribed operating domain. 𝐹𝐴, 𝐹𝐵 , 𝐹2𝑙 , and 𝐹1𝑙 were sampled within ±15% of their nominal values, while 𝑇𝐴 and 𝑇𝐵 were sampled over 305–325 K and 290–310 K, respectively. For each sampled input combination, the DAE was solved for steady state conditions using the IDA DAE solver from scikit-SUNDAE. Parameter continuation was employed by moving from the nominal operating point to each sampled condition through 12 intermediate steps, with the solution from each step used as the initial guess for the subsequent step. A solution was retained only when the IDA calculation converged, all predicted variables were finite, and the maximum absolute residual of the 19 DAE equations was less than 10−7 . To keep the problem underdefined, we do not consider the steady state conditions as constraints in the projection layer. The following thermodynamic constraints are considered as equality constraints: 0 = 𝑁 𝐴1 − 𝑘 𝐴𝑎(𝑦 𝐴1 − 𝑦 ∗𝐴1 )

𝑀1𝑙

𝑀1𝑙 , 𝜌 𝑀𝑙 ∗ 0 = 𝑁𝐶2 − 𝑘 𝐶 𝑎(1 − 𝑦 𝐴2 − 𝑦 𝐶2 ) 2, 𝜌 ! 𝑙 𝑀 0 = 𝑃2 𝑉2𝑇 − 2 − 𝑀2𝑣 𝑅𝑇2 , 𝜌 ∗ 0 = 𝑁𝐶1 − 𝑘 𝐶 𝑎 𝑦 𝐶1 − (1 − 𝑦 𝐴1 )

,

𝜌 𝑀2𝑙

0 = 𝑁 𝐴2 − 𝑘 𝐴𝑎(𝑦 ∗𝐴2 − 𝑦 ∗𝐴2 ) , 𝜌 ! 𝑀1𝑙 − 𝑀1𝑣 𝑅𝑇1 , 0 = 𝑃1 𝑉1𝑇 − 𝜌 0 = 𝑃1 − 𝑃2 −

1 (𝐹1𝑣 ) 7/4 . 0.09

(7a) (7b) (7c) (7d)

13

Only interphase mass transfer and pressure drop correlations are considered in the projection layer to keep the system underdefined. Additionally, the following inequality constraints are applied: 𝑀1𝑙 ∈ [5.0, 44.9]

𝑀1𝑣 ∈ [0.01, 5.0]

𝑀2𝑙 ∈ [0.1, 44.0] 𝑀2𝑣 ∈ [0.3, 6.5] 𝑥 𝐴1 ∈ [0.35, 0.70] 𝑥 𝐵1 ∈ [0.18, 0.55] 𝑥 𝐴1 + 𝑥 𝐵1 <= 1.0 𝑥 𝐴2 ∈ [0.55, 0.82] 𝑦 𝐴1 ∈ [0.35, 0.56] 𝑦 𝐴2 ∈ [0.22, 0.46] 𝑇1 ∈ [315.0, 345.0] 𝑇2 ∈ [280.0, 315.0] 𝑃1 ∈ [45.0, 55.0] 𝑃2 ∈ [45.0, 55.0] 𝐹1𝑣 ∈ [0.15, 0.27] 𝐹2𝑣 ∈ [0.10, 0.18] 𝑁 𝐴1 ∈ [−0.02, 0.02] 𝑁𝐶1 ∈ [0.09, 0.13] 𝑁 𝐴2 ∈ [0.03, 0.07] 𝑁𝐶2 ∈ [0.01, 0.03]

(8)

to ensure an interior point as initialization from the backbone. We observe that without these inequality constraints the backbone produces a highly infeasible initialization which often diverges from the desired solution during projection. We use 1000 different parameter conditions and split the dataset 85 : 15 ratio for training and validation. The network backbone structure is given as 128 hidden nodes and 2 hidden layers with ReLU activation. We use “Fischer-Burmeister-Reduced" formulation of the distance minimization problem with tolerance 10−9 , maximum iterations 50 and smoothing 10−9 along with other default configurations given in the code. We train the model for 3000 epochs with a batch size of 50, and learning rate 10−3 . For the first 2000 epochs we pretrain the backbone network for a suitable initialization to the projection. The training is supervised where we fit the data generated in the procedure mentioned before. We consider 7 equality constraints 7 and 20 inequalities 8 in the projection layer. The code to train the model is given in codeblock 2. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32

import os from pathlib import Path import jax import jax.numpy as jnp os.environ["JAX_ENABLE_X64"] = "True" jax.config.update("jax_enable_x64", True) from kkthn import KKTHardNet, Loss, MLP, Projection SEED = 42 # Configurations TRAIN = {"epochs": 700, "epoch_mlp": 600, "batch_size": 50, "learning_rate": 1e-3, "train_frac": 0.85, "seed": SEED, "dtype": "float64", "print_every": 100, "minimum_projection_convergence": 0.95, "readiness_extra_epochs": 500, "readiness_check_every": 10, "backend": "cpu", "output_dir": str(Path.cwd() / "runs")} # Parameters and Variables PARAMETER_NAMES = ["FA", "FB", "TA", "TB", "F2l", "F1l"] VARIABLE_NAMES = ["M1l", "xA1", "xB1", "M1v", "yA1", "T1", "M2l", "xA2", "M2v", "yA2", "T2", "NA1", "NC1", "NA2", "NC2", "P1", "P2", "F1v", "F2v"] DATA_PATH = Path.cwd() / "reactor_condensor" / "reactor_condensor.csv" # Backbone nn = MLP(hidden_sizes=(128, 128), activation="tanh", predict_duals=True) # Projection Layer projection = Projection(name="reactor_condensor", method="fischer-burmeister-reduced", config={"smoothing": 1e-9, "tolerance": 1e-9, "max_iterations": 50}) x = projection.add_parameter(PARAMETER_NAMES)

14

33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91

y = projection.add_variable(VARIABLE_NAMES) a, kA, kC, rho = 1000.0, 20.0, 30.0, 15.0 V1T, V2T, R_gas = 3.0, 3.0, 0.082057 MOL_TO_KMOL, PA_PER_ATM = 1.0e-3, 101325.0 PsA1 = projection.exp(25.1 - 3400.0 / (y.T1 + 20.0)) / PA_PER_ATM PsC1 = projection.exp(27.3 - 4100.0 / (y.T1 + 70.0)) / PA_PER_ATM PsA2 = projection.exp(25.1 - 3400.0 / (y.T2 + 20.0)) / PA_PER_ATM PsC2 = projection.exp(27.3 - 4100.0 / (y.T2 + 70.0)) / PA_PER_ATM yA1_star = PsA1 * y.xA1 / y.P1 yC1_star = PsC1 * (1.0 - y.xA1 - y.xB1) / y.P1 yA2_star = PsA2 * y.xA2 / y.P2 yC2_star = PsC2 * (1.0 - y.xA2) / y.P2 projection.constraints.add( y.NA1 - kA * a * (y.yA1 - yA1_star) * (y.M1l / rho) * MOL_TO_KMOL == 0, y.NC1 - kC * a * (yC1_star - (1.0 - y.yA1)) * (y.M1l / rho) * MOL_TO_KMOL == 0, y.NA2 - kA * a * (y.yA2 - yA2_star) * (y.M2l / rho) * MOL_TO_KMOL == 0, y.NC2 - kC * a * (1.0 - y.yA2 - yC2_star) * (y.M2l / rho) * MOL_TO_KMOL == 0, y.P1 * (V1T - y.M1l / rho) - y.M1v * R_gas * y.T1 == 0, y.P2 * (V2T - y.M2l / rho) - y.M2v * R_gas * y.T2 == 0, y.P1 - y.P2 - (1.0 / 0.09) * y.F1v ** (7.0 / 4.0) == 0 ) projection.constraints.add( y.M1l >= 5.0, y.M1l <= 44.9, y.M1v >= 0.01, y.M1v <= 5.0, y.M2l >= 0.1, y.M2l <= 44.0, y.M2v >= 0.3, y.M2v <= 6.5, y.xA1 >= 0.35, y.xA1 <= 0.70, y.xB1 >= 0.18, y.xB1 <= 0.55, y.xA1 + y.xB1 <= 1.0, y.xA2 >= 0.55, y.xA2 <= 0.82, y.yA1 >= 0.35, y.yA1 <= 0.56, y.yA2 >= 0.22, y.yA2 <= 0.46, y.T1 >= 315.0, y.T1 <= 345.0, y.T2 >= 280.0, y.T2 <= 315.0, y.P1 >= 45.0, y.P1 <= 55.0, y.P2 >= 45.0, y.P2 <= 55.0, y.F1v >= 0.15, y.F1v <= 0.27, y.F2v >= 0.10, y.F2v <= 0.18, y.NA1 >= -0.02, y.NA1 <= 0.02, y.NC1 >= 0.09, y.NC1 <= 0.13, y.NA2 >= 0.03, y.NA2 <= 0.07, y.NC2 >= 0.01, y.NC2 <= 0.03 ) # Loss function loss = Loss.mse(primal_consistency_weight=1.0, equality_dual_consistency_weight=0.0, inequality_dual_consistency_weight=0.0) # Model definition model = KKTHardNet("surrogate", nn, loss, projection, train=TRAIN) # Dataset model.dataset(parameters=DATA_PATH, variables=DATA_PATH) # Training fb_reduced_result = model.model() model.summary() model.plot_history()

# Inference model = KKTHardNet().load(Path("/path/to/metadata.json")) x_new = jnp.array([[0.08981, 0.0533, 305.206, 308.4160856, 0.0678, 0.0108]]) y_pred = model.predict(x_new)

Code 2: Code for training the Reactor-Condensor model with KKT-HardNet.

15

Table 5: Performance comparison for the reactor–condenser problem. We define the metrics in appendix A. MSE MLP PINN ENFORCE KKT-HardNet

9.75 × 10−2 1.31 × 10−1 7.21 × 100 3.11 × 10−3

Mean eq. 0.29 0.61 7.36 × 10−10 1.33 × 10−7

Max eq.

Mean ineq.

Max ineq.

Training Time (s)

Inference Time (ms/batch)

7.429 5.29 5.10 × 10−9 8.59 × 10−7

2.22 × 10−5

6.73 × 10−2

7.03 × 10−3 0.00 6.19 × 10−8

9.80 × 10−1 0.00 9.45 × 10−7

39.78 41.00 16362.92 786.21

0.35 0.31 208.39 57.86

Table 5 compares the surrogate accuracy, constraint satisfaction, and computational performance of the different learning approaches. KKT-HardNet achieves the lowest prediction error, with an MSE of 3.11×10−3 , compared to other approached. At the same time, the projected predictions maintain small constraint residuals, with maximum equality and inequality violations of 8.59 × 10−7 and 9.45 × 10−7 , respectively. Although ENFORCE attains tighter constraint residuals, it incurs substantially higher training and inference costs, requiring 16362.92 s for training and 208.39 ms per inference batch, compared to KKT-HardNet. KKTHardNet takes only 786.21 s and 57.86 ms for training and inference, respectively. The unconstrained MLP and soft-constrained PINN are computationally cheaper, but exhibit considerably larger physical-constraint violations, demonstrating the trade-off between computational cost, predictive accuracy, and strict physical consistency.

4

PCML for Real-Time Optimization

Real-time optimization (RTO) determines economically optimal process decisions in response to changing process and economic conditions. Quantities such as energy prices, feed conditions, production targets, and equipment availability can vary with time, requiring the preferred operating point to change. RTO addresses this problem by repeatedly solving a constrained optimization problem while accounting for process physics and operating limits. In general, let 𝜽 ∈ R𝑛 𝜃 denote the time-varying process or economic parameters and 𝒙 ∈ R𝑛𝑥 the decision variables. The RTO problem can be expressed as 𝒙 ∗ (𝜽) = arg min 𝒙 𝑓 (𝜽, 𝒙) subject to ℎ(𝜽, 𝒙) = 0 and 𝑔(𝜽, 𝒙) ≤ 0, where 𝑓 is the economic objective and ℎ and 𝑔 represent the process and operating constraints. Repeatedly solving such optimization problems can become computationally expensive, particularly for nonlinear systems and applications requiring frequent updates. KKT-HardNet provides an alternative by learning the parametric solution map 𝜽 ↦→ 𝒙 ∗ offline. During training, the loss function drives the network toward economically favorable solutions, while the projection layer enforces the constraints. After training, changes in the process or economic parameters can be mapped directly to feasible near-optimal operating conditions through model inference, reducing the computational burden associated with repeatedly solving the original optimization problem online. The general RTO formulation with KKT-HardNet is illustrated in Figure 9. 4.1

A CSTR Optimization Case Study

We consider a continuous stirred tank reactor (CSTR) with the reversible exothermic reaction 𝐴 ↔ 𝐵 [39]. The process states are the concentrations 𝐶 𝐴 and 𝐶 𝐵 and the reactor temperature 𝑇 , while 𝑄 denotes the heat supplied through the reactor jacket. The nonlinear mass and energy balances are as follows: 𝑑𝐶 𝐴 1 = (𝐶 𝐴0 − 𝐶 𝐴) − 𝑘 𝐴 𝑒 −𝐸 𝐴/(𝑅𝑇 ) 𝐶 𝐴 + 𝑘 𝐵 𝑒 −𝐸𝐵 /(𝑅𝑇 ) 𝐶 𝐵 , 𝑑𝑡 𝜏 𝑑𝐶 𝐵 1 = − 𝐶 𝐵 + 𝑘 𝐴 𝑒 −𝐸 𝐴/(𝑅𝑇 ) 𝐶 𝐴 − 𝑘 𝐵 𝑒 −𝐸𝐵 /(𝑅𝑇 ) 𝐶 𝐵 , 𝑑𝑡 𝜏  1 𝑑𝑇 Δ𝐻  𝑄 =− 𝑘 𝐴 𝑒 −𝐸 𝐴/(𝑅𝑇 ) 𝐶 𝐴 − 𝑘 𝐵 𝑒 −𝐸𝐵 /(𝑅𝑇 ) 𝐶 𝐵 + (𝑇0 − 𝑇) + . 𝑑𝑡 𝜌𝐶 𝑃 𝜏 𝜌𝐶 𝑃 𝑉

16

(9)

min x,𝜃

f(𝒙; 𝜽)

s. t

𝒉(𝒙; 𝜽) = 0 𝒈(𝒙; 𝜽) ≤ 0 𝐿obj

Inputs 𝜽

..

𝐿𝑐𝑜𝑛𝑠𝑖𝑠𝑡𝑒𝑛𝑐𝑦

𝒙ƶ

Outputs 𝒙ƿ

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 𝒙˜ .

Here, 𝑇0 and 𝐶 𝐴0 denote the feed temperature and concentration; 𝑘 𝐴 and 𝑘 𝐵 are the pre-exponential factors; 𝐸 𝐴 and 𝐸 𝐵 are the activation energies; 𝜏 is the residence time; and Δ𝐻 , 𝜌 , 𝐶 𝑃 , and 𝑉 denote the reaction enthalpy, liquid density, heat capacity, and reactor volume, respectively. The model parameters are summarized in Table 6.

Table 6: Parameter values and steady-state values for the CSTR case study. 𝑇0 = 400 K 𝑘 𝐴 = 5000 /s 𝐸 𝐴 = 1 × 104 cal/mol 𝑅 = 1.987 cal/(mol K) 𝜌 = 1 kg/L 𝐶 𝐴0 = 1 mol/L 𝐶 𝐴𝑠 = 426.743 K

𝜏 = 60 s 𝑘 𝐵 = 106 /s 𝐸 𝐵 = 1.5 × 104 cal/mol Δ𝐻 = −5000 cal/mol 𝐶 𝑃 = 1000 cal/(kg K) 𝑉 = 100 L 𝑄 𝑠 = 40386 cal/s

Because the reaction is reversible, excessively low temperature results in slow conversion, whereas high temperature increasingly favors the reverse reaction. The economically preferred operating point, therefore, depends on both reactant conversion and the cost of heat. We define the steady-state economic objective as

𝐽RTO (𝐶 𝐴, 𝑄; 𝑝) =

𝐶𝐴 + 𝑝𝑄, 𝐶 𝐴0

where 𝑝 is the time-varying heat price. The corresponding parametric constrained RTO problem is 17

(10)

min 𝐶 𝐴 ,𝐶𝐵 ,𝑇 ,𝑄

s.t.

𝐶𝐴 + 𝑝𝑄 𝐶 𝐴0 1 (𝐶 𝐴0 − 𝐶 𝐴) − 𝑘 𝐴 𝑒 −𝐸 𝐴/(𝑅𝑇 ) 𝐶 𝐴 + 𝑘 𝐵 𝑒 −𝐸𝐵 /(𝑅𝑇 ) 𝐶 𝐵 = 0, 𝜏 1 − 𝐶 𝐵 + 𝑘 𝐴 𝑒 −𝐸 𝐴/(𝑅𝑇 ) 𝐶 𝐴 − 𝑘 𝐵 𝑒 −𝐸𝐵 /(𝑅𝑇 ) 𝐶 𝐵 = 0, 𝜏  1 Δ𝐻  𝑄 − 𝑘 𝐴 𝑒 −𝐸 𝐴/(𝑅𝑇 ) 𝐶 𝐴 − 𝑘 𝐵 𝑒 −𝐸𝐵 /(𝑅𝑇 ) 𝐶 𝐵 + (𝑇0 − 𝑇) + = 0, 𝜌𝐶 𝑃 𝜏 𝜌𝐶 𝑃 𝑉 0 ≤ 𝐶 𝐴 ≤ 1, 0 ≤ 𝐶 𝐵 ≤ 1, 400 ≤ 𝑇 ≤ 500, 0 ≤ 𝑄 ≤ 105 .

(11)

KKT-HardNet is trained directly as an optimization proxy for Eq. 11, with heat price 𝑝 as the problem parameter (𝜃 ) and (𝐶 𝐴, 𝐶 𝐵 , 𝑇, 𝑄) as the optimization variables (𝑥 ). We sample 2500 heat prices uniformly over 𝑝 ∈ [10−7 , 2 × 10−6 ] and use an 80:20 training–validation split. The backbone consists of three hidden layers with 32 neurons per layer and ReLU activation. The Fischer–Burmeister-Reduced projection is used with tolerance 10−10 , maximum 300 iterations, and smoothing parameter 10−10 . The model is trained for 2000 epochs using a batch size of 50, learning rate 3 × 10−4 , and primal consistency weight (𝜆 𝑐 ) to be 0.05. The loss function is unsupervised and no label data is required. The implementation is shown in Code 3. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40

import os os.environ["JAX_ENABLE_X64"] = "True" from pathlib import Path import jax jax.config.update("jax_enable_x64", True) import jax.numpy as jnp from kkthn import KKTHardNet, Loss, MLP, Projection # System Parameters SEED = 50, T_FEED = 400.0, CA_FEED = 1.0, TAU = 60.0, KA = 5000.0, KB = 1.0e6, EA = 1.0e4, EB = 1.5e4, R = 1.987, DEL_H = -5000.0, RHO = 1.0, CP = 1000.0, V = 100.0 # Constructor blocks def equality_constraints(projection, y): rA = KA * projection.exp(-EA/(R*y.T_s))*y.CA_s rB = KB * projection.exp(-EB/(R*y.T_s))*y.CB_s r = rA-rB return [(CA_FEED-y.CA_s)/TAU-r==0.0, -y.CB_s/TAU+r==0.0, (-DEL_H/(RHO*CP)) * r + (T_FEED-y.T_s) / TAU + y.Q_s/(RHO*CP*V) == 0.0] def inequality_constraints(y): return [y.CA_s >= 0.0, y.CA_s <= 1.0, y.CB_s >= 0.0, y.CB_s <= 1.0, y.T_s >= 400.0, y.T_s <= 500.0, y.Q_s >= 0.0, y.Q_s <= 1.0e5] def objective_function(x, y): return y.CA_s / CA_FEED + x.heat_price * y.Q_s # Configuration TRAIN = { "epochs": 2000, "batch_size": 50, "learning_rate": 3e-4, "train_frac": 0.8, "seed": SEED, "dtype": "float64", "print_every": 10, "standardize_inputs": True, "pretrain": False, "drop_last_batch": True, "backend": "cpu", "output_dir": str(Path.cwd() / "runs")} DATA_PATH = Path.cwd() / "rto" / "heat_price.csv" # Backbone nn = MLP(hidden_sizes=(32, 32, 32), activation="relu", predict_duals=True, inequality_dual_floor=1e-8) # Projection Layer projection = Projection(name="CSTR_RTO", method="fischer-burmeister-reduced",

18

41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

config={"tolerance": 1e-10, "max_iterations": 300, "smoothing": 1e-10}) x = projection.add_parameter(["heat_price"]) y = projection.add_variable(["CA_s", "CB_s", "T_s", "Q_s"]) projection.constraints.add(equality_constraints(projection, y)) projection.constraints.add(inequality_constraints(y)) projection.objective = objective_function(x, y) # Loss function loss = Loss.objective(primal_consistency_weight=0.05, equality_dual_consistency_weight=0.00, inequality_dual_consistency_weight=0.00) # Model and Training model = KKTHardNet("optimization", nn, loss, projection, train=TRAIN) model.dataset(parameters=DATA_PATH) fb_reduced_result = model.optimize() # Inference model = KKTHardNet().load(Path("/path/to/metadata.json")) x_new = jnp.array([[1.0e-6]]) y_pred = model.predict(x_new)

Code 3: Code for training the real-time optimization (RTO) model with KKT-HardNet.

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. IPOPT PINN DC3 KKT-HardNet

Objective

Mean eq.

Max eq.

Mean ineq.

Max ineq.

Training Time (s)

Inference Time (ms/batch)

0.535 −0.487 0.560 0.535

3.43 × 10−12

9.98 × 10−11

1.96 × 10−2

4.06 × 10−2

1.11 × 10−16 4.04 × 10−8

5.71 × 10−18 8.57 × 10−8

0.00 0.06 0.00 0.00

0.00 0.514 0.00 0.00

−− 56.87 333.38 3725.66

1.83∗ 0.086 2.13 1.65

*Solver time is reported as average solve time considering full parallelization across 100 instances

As shown in Table 7, KKT-HardNet reproduces the reference optimizer objective value of 0.535 while maintaining very small tolerance for equality residuals and zero inequality violation. The PINN obtains an apparently lower objective but violates both equality and inequality constraints, whereas DC3 remains feasible but does not find the optimal objective value. KKT-HardNet also achieves smaller inference time of 1.65 ms for a batch of 100 instances, compared to 1.83 ms for the reference optimizer. The learned RTO policy is further evaluated over the complete heat-price range, as shown in Fig. 10. Increasing the heat price reduces the optimal heating rate from approximately 5.62 × 104 to 1.33 × 104 cal/s and decreases the optimal reactor temperature from 436.27 to 410.34 K. Correspondingly, 𝐶 𝐴𝑠 𝑝 increases from 0.491 to 0.533 mol/L and 𝐶 𝐵𝑠 𝑝 decreases from 0.509 to 0.467 mol/L, reflecting the trade-off between conversion and energy cost. The optimal objective increases from 0.497 to 0.560. The sensitivity profiles show the same trend, with d𝑄 𝑠 𝑝 /d𝑝 < 0, d𝑇 𝑠 𝑝 /d𝑝 < 0, d𝐶 𝐴𝑠 𝑝 /d𝑝 > 0, and d𝐶 𝐵𝑠 𝑝 /d𝑝 < 0. For real-time evaluation, the trained RTO model was tested over 500 warm single-instance predictions. The mean and median inference latencies were 38.34 and 38.32 ms, respectively, with a maximum of 66.49 ms, which is negligible compared to the 1000 s RTO update interval. The first model evaluation required 7.67 s because of one-time JAX compilation and was excluded from the online timing statistics.

5

PCML-based Model Predictive Control

Model predictive control (MPC) determines control actions by repeatedly solving a finite-horizon optimization problem using the current process state while explicitly accounting for system dynamics and operating 19

(a) Economic operating map over the complete heat-price range.

(b) Sensitivity of the optimal operating point with respect to heat price.

Figure 10: Learned RTO policy and its local sensitivity to heat price. constraints. At each sampling instance, a sequence of future states and control actions is optimized, the first control action is applied to the process, and the optimization is repeated using the updated state. In general, let 𝒙 𝑘 denote the current state and other time-varying controller parameters, and 𝒚 contain the predicted states and control actions over the horizon. The MPC problem can be expressed as 𝒚 ∗ (𝒙 𝑘 ) = arg min𝒚 𝐽 (𝒙 𝑘 , 𝒚) subject to ℎ(𝒙 𝑘 , 𝒚) = 0 and 𝑔(𝒙 𝑘 , 𝒚) ≤ 0, where 𝐽 is the finite-horizon control objective and ℎ and 𝑔 represent the discrete time governing equations and operating constraints. Since the optimization must be solved at every sampling instance, the computational cost of nonlinear MPC can become significant for systems with short control intervals or large prediction horizons. KKT-HardNet can instead learn the parametric solution map 𝒙 𝑘 ↦→ 𝒚 ∗ offline, while incorporating the MPC constraints through the projection layer. During online operation, the current state and controller parameters are supplied to the trained model to obtain a feasible near-optimal control solution through inference, from which only the first control action is implemented according to the receding-horizon principle. This provides a general approach for accelerating repeated constrained optimization in MPC while retaining the specified process and operating constraints. The general formulation of real-time optimization with KKT-HardNet is illustrated in Figure 11. 5.1

A CSTR Control Example

We consider the same CSTR system for the control case study. Given the optimal steady-state condition (𝑥 𝑠 𝑝 , 𝑄 𝑠 𝑝 ) , with 𝑥 𝑠 𝑝 = [𝐶 𝐴𝑠 𝑝 , 𝐶 𝐵𝑠 𝑝 , 𝑇 𝑠 𝑝 ] ⊤ , the MPC uses the current state 𝑥 = [𝐶 𝐴, 𝐶 𝐵 , 𝑇] ⊤ together with these set points to determine the heat input 𝑄 . The MPC is formulated in deviation coordinates, 𝑥¯ = 𝑥 − 𝑥 𝑠 𝑝 and 𝑢 = 𝑄 − 𝑄 𝑠 𝑝 , such that 𝑥 = 𝑥 𝑠 𝑝 + 𝑥¯ and 𝑄 = 𝑄 𝑠 𝑝 + 𝑢 . Defining the net reaction rate as 𝑟 = 𝑘 𝐴 𝑒 −𝐸 𝐴/(𝑅𝑇 ) 𝐶 𝐴 − 𝑘 𝐵 𝑒 −𝐸𝐵 /(𝑅𝑇 ) 𝐶 𝐵 , the CSTR dynamics can be written in control-affine form as 𝑥¤ = 𝑓 (𝑥; 𝑄 𝑠 𝑝 ) + 𝑔𝑢,

where 20

(12)

𝐿obj

Inputs 𝒙(𝒕𝒌 ), 𝒙𝒔𝒑 , 𝜽

𝐿𝑐𝑜𝑛𝑠𝑖𝑠𝑡𝑒𝑛𝑐𝑦

Outputs ഥ𝑻𝟏|𝒌 , . . . , 𝒙 ഥ𝑻𝑵|𝒌 𝒙 ഥ 𝑻𝟎|𝒌 , . . . , 𝒖 ഥ 𝑻𝑵−𝟏|𝒌 𝒖 k = 1, 2, . . . , N

..

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. 𝐶 𝐴0 − 𝐶 𝐴     −𝑟   𝜏   𝐶𝐵 𝑠𝑝 ,  − +𝑟 𝑓 (𝑥; 𝑄 ) =   𝜏  Δ𝐻  𝑠 𝑝 𝑇0 − 𝑇 𝑄   𝑟 + + −   𝜏 𝜌𝐶 𝑃 𝑉   𝜌𝐶 𝑃

 0     0  𝑔=  1 .      𝜌𝐶 𝑃 𝑉 

(13)

To promote stabilizing closed-loop behavior, we consider the Lyapunov function 𝑉 ( 𝑥) ¯ = 105𝐶¯ 2𝐴 +105𝐶¯ 𝐵2 + 𝑇¯ 2 , ⊤ ⊤ with 𝐿 𝑓 𝑉 = ∇𝑉 𝑓 and 𝐿 𝑔𝑉 = ∇𝑉 𝑔. The auxiliary stabilizing control law is 

𝐿 𝑓𝑉 +

√︁

ℎ( 𝑥) ¯ =−

 (𝐿 𝑓 𝑉) 2 + (𝐿 𝑔𝑉) 4 𝐿 𝑔𝑉 (𝐿 𝑔𝑉) 2 + 𝜖

,

𝜖 = 10−12 .

(14)

Defining 𝑉¤ ( 𝑥, ¯ 𝑢) = 𝐿 𝑓 𝑉 ( 𝑥) ¯ + 𝐿 𝑔𝑉 ( 𝑥)𝑢 ¯ , the first control action is constrained such that 𝑉¤ ( 𝑥¯0| 𝑘 , 𝑢 0| 𝑘 ) ≤ ¤ 𝑉 ( 𝑥¯0| 𝑘 , ℎ( 𝑥¯0| 𝑘 )) . The state bounds are 𝑥 min = [0, 0, 400] ⊤ and 𝑥 max = [1, 1, 500] ⊤ , the heat-input bounds are 𝑄 min = 0 and 𝑄 max = 105 , and consecutive predicted control moves are limited by Δ𝑄 max = 5 × 103 . The resulting parametric MPC problem is as follows: min

𝑥¯ 1:𝑁 |𝑘 , 𝑢0:𝑁 −1|𝑘

s.t.

Δ𝑡

𝑁 ∑︁

∥ 𝑥¯ 𝑗 | 𝑘 ∥ 2𝑄 𝑥 + Δ𝑡

𝑗=1

𝑁 −1 ∑︁

𝑅𝑢 𝑢 2𝑗 | 𝑘

𝑗=0 𝑠𝑝

𝑥¯0| 𝑘 = 𝑥(𝑡 𝑘 ) − 𝑥 ,

  𝑥¯ 𝑗+1| 𝑘 = 𝑥¯ 𝑗 | 𝑘 + Δ𝑡 𝑓 (𝑥 𝑠 𝑝 + 𝑥¯ 𝑗 | 𝑘 ; 𝑄 𝑠 𝑝 ) + 𝑔𝑢 𝑗 | 𝑘

𝑥 min ≤ 𝑥 𝑠 𝑝 + 𝑥¯ 𝑗 | 𝑘 ≤ 𝑥max , 𝑄 min ≤ 𝑄 𝑠 𝑝 + 𝑢 𝑗 | 𝑘 ≤ 𝑄 max , 𝑉¤ ( 𝑥¯0| 𝑘 , 𝑢 0| 𝑘 ) ≤ 𝑉¤ ( 𝑥¯0| 𝑘 , ℎ( 𝑥¯0| 𝑘 )).

(15)

|𝑢 𝑗 | 𝑘 − 𝑢 𝑗 −1| 𝑘 | ≤ Δ𝑄 max

The problem is parameterized by 𝜃 𝑘 = [𝐶 𝐴 (𝑡 𝑘 ), 𝐶 𝐵 (𝑡 𝑘 ), 𝑇 (𝑡 𝑘 ), 𝐶 𝐴𝑠 𝑝 , 𝐶 𝐵𝑠 𝑝 , 𝑇 𝑠 𝑝 , 𝑄 𝑠 𝑝 ] ⊤ . Only the first optimized control action is implemented, giving 𝑄(𝑡 𝑘 ) = 𝑄 𝑠 𝑝 + 𝑢 ∗0| 𝑘 , after which the process is propagated to the next sampling instance and the optimization is repeated. 21

We use a sampling time of 5 s and a prediction/control horizon of 50 s, corresponding to 𝑁 = 10. The resulting NLP contains 7 parameters, 40 primal variables, 30 equality constraints, and 99 inequality constraints. The 2500 RTO operating points are used to construct the MPC training conditions with an 80:20 training– validation split. The neural-network backbone contains three hidden layers with 256 neurons per layer and ReLU activation. The projection uses the primal-dual interior-point method with tolerance 10−8 , maximum 300 iterations, and barrier regularization parameter 𝜇 = 10−4 . The model is trained for 1200 epochs using a batch size of 50, learning rate 3 × 10−4 , and primal consistency weight 0.01. The implementation is shown in Code 4. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52

import os os.environ["JAX_ENABLE_X64"] = "True" from pathlib import Path import jax jax.config.update("jax_enable_x64", True) import jax.numpy as jnp from kkthn import KKTHardNet, Expression, Loss, MLP, Projection

SEED = 42, T_FEED = 400.0, CA_FEED = 1.0, TAU = 60.0, KA = 5000.0, KB = 1.0e6, EA = 1.0e4, EB = 1.5e4, R = 1.987, DEL_H = -5000.0, RHO = 1.0, CP = 1000.0, V = 100.0 DT = 5.0, N = 10, QC_CA = 1.0, QC_CB = 1.0, QC_T = 5e-5, RC = 1e-11, Q_MIN = 0.0, Q_MAX = 1.0e5, DU_MAX = 5.0e3 PARAMETER_NAMES = ["C_A0", "C_B0", "T0", "C_Asp", "C_Bsp", "Tsp", "Qsp"] STATE_VARIABLE_NAMES = [f"x{v}_{k}" for k in range(1, N+1) for v in ("CA", "CB", "T")] CONTROL_VARIABLE_NAMES = [f"u{k}" for k in range(N)] VARIABLE_NAMES=STATE_VARIABLE_NAMES + CONTROL_VARIABLE_NAMES def objective(theta, y): return DT * sum(QC_CA*y[3*k]**2 + QC_CB*y[3*k+1]**2 + QC_T*y[3*k+2]**2 + RC*y[30+k]**2 for k in range(N)) def equalities(theta, y): CA0, CB0, T0, CAsp, CBsp, Tsp, Qsp = theta residuals = [] for k in range(N): prev = ((CA0 - CAsp, CB0 - CBsp, T0 - Tsp) if k == 0 else y[3 * (k - 1):3 * k]) xCA_prev, xCB_prev, xT_prev = prev xCA_next, xCB_next, xT_next = y[3 * k:3 * k + 3] u = y[30 + k] CA = CAsp + xCA_prev CB = CBsp + xCB_prev T = Tsp + xT_prev Q = Qsp + u r = (KA * jnp.exp(-EA / (R * T)) * CA - KB * jnp.exp(-EB / (R * T)) * CB) residuals += [xCA_next - xCA_prev - DT * ((CA_FEED - CA) / TAU - r), xCB_next - xCB_prev - DT * (-CB / TAU + r), xT_next - xT_prev - DT*((-DEL_H/(RHO*CP))*r + (T_FEED-T)/TAU + Q/(RHO*CP*V))] return jnp.stack(residuals) def inequalities(theta, y): CA0, CB0, T0, CAsp, CBsp, Tsp, Qsp = theta xCA0, xCB0, xT0 = CA0 - CAsp, CB0 - CBsp, T0 - Tsp rA0 = KA * jnp.exp(-EA / (R * T0)) * CA0 rB0 = KB * jnp.exp(-EB / (R * T0)) * CB0 f1 = (CA_FEED - CA0) / TAU - rA0 + rB0 f2 = -CB0 / TAU + rA0 - rB0 f3 = (-DEL_H / (RHO * CP)) * (rA0 - rB0) + (T_FEED - T0) / TAU + Qsp / (RHO * CP * V) V_CA, V_CB, V_T = 2.0e5 * xCA0, 2.0e5 * xCB0, 2.0 * xT0 LfV = V_CA * f1 + V_CB * f2 + V_T * f3 LgV = V_T / (RHO * CP * V) h = -(LfV + jnp.sqrt(LfV**2 + LgV**4)) * LgV / (LgV**2 + 1e-12)

22

53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116

residuals = [LfV + LgV * y[30] - (LfV + LgV * h)] for k in range(N): xCA, xCB, xT = y[3*k:3*k+3] u = y[30+k] residuals += [-xCA-CAsp, xCA+CAsp-1, -xCB-CBsp, xCB+CBsp-1, 400-Tsp-xT, xT+Tsp-500, Q_MIN-Qsp-u, u+Qsp-Q_MAX] if k: du = u-y[30+k-1] residuals += [du-DU_MAX, -du-DU_MAX] return jnp.stack(residuals) def _theta_from_ctx(ctx): return jnp.stack([ctx.parameters[ctx.parameter_index[n]] for n in PARAMETER_NAMES]) def lmpc_primal_transform(raw_y, X): temperature_indices = jnp.asarray([2, 5, 8, 11, 14, 17, 20, 23, 26, 29]) Tsp = X[..., 5:6] T = Tsp + raw_y[..., temperature_indices] xT = 400.0 + jax.nn.relu(T - 400.0) - Tsp return raw_y.at[..., temperature_indices].set(xT) # Configurations TRAIN = {"epochs": 1200, "batch_size": 50, "learning_rate": 3e-4, "train_frac": 0.8, "seed": SEED, "dtype": "float64", "print_every": 1, "standardize_inputs": True, "pretrain": False, "drop_last_batch": True, "backend": "cpu", "output_dir": str(Path.cwd()/"runs")} DATA_PATH = Path.cwd()/"lmpc"/"parameters_new.csv" # Backbone nn = MLP(hidden_sizes=(256, 256, 256), activation="relu", predict_duals=True, inequality_dual_floor=1e-8) # Projection Layer projection = Projection(name="CSTR_LMPC", method="primal-dual-interior-point", config={"tolerance": 1e-8,"max_iterations": 300,"mu": 1e-4}) projection.add_parameter(PARAMETER_NAMES) projection.add_variable(VARIABLE_NAMES) equality_dependencies, inequality_dependencies = [], [("u0",)] for k in range(N): previous = [] if k == 0 else [f"xCA_{k}", f"xCB_{k}", f"xT_{k}"] for v in ("xCA", "xCB", "xT"): equality_dependencies.append( tuple(previous + [f"{v}_{k+1}"] + ([f"u{k}"] if v == "xT" else []))) for v in ("xCA", "xCB", "xT"): inequality_dependencies += [(f"{v}_{k+1}",), (f"{v}_{k+1}",)] inequality_dependencies += [(f"u{k}",), (f"u{k}",)] if k > 0: inequality_dependencies += [(f"u{k-1}", f"u{k}"),(f"u{k-1}", f"u{k}")] for i, deps in enumerate(equality_dependencies): expr = Expression( lambda ctx, i=i:equalities(_theta_from_ctx(ctx), ctx.variables)[i],f"h[{i}]", variable_dependencies=deps,parameter_dependencies=PARAMETER_NAMES, affine_in_variables=False,) projection.constraints.equality(expr) for i, deps in enumerate(inequality_dependencies): expr = Expression( lambda ctx, i=i:inequalities(_theta_from_ctx(ctx), ctx.variables)[i],f"g[{i}]", variable_dependencies=deps,parameter_dependencies=PARAMETER_NAMES, affine_in_variables=True,) projection.constraints.inequality(expr) projection.objective = Expression(

23

117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135

lambda ctx: objective(_theta_from_ctx(ctx),ctx.variables,),"objective", variable_dependencies=VARIABLE_NAMES, parameter_dependencies=PARAMETER_NAMES, affine_in_variables=False,) # Loss function loss = Loss.objective(primal_consistency_weight=0.01, equality_dual_consistency_weight=0.0, inequality_dual_consistency_weight=0.0) # Model definition model = KKTHardNet("optimization", nn, loss, projection, train=TRAIN) model.primal_transform_fn = lmpc_primal_transform model.dataset(parameters=DATA_PATH) # Training pdip_result = model.optimize() # Inference model = KKTHardNet().load(Path("/path/to/metadata.json")) x_new = jnp.array([[0.50, 0.45, 425.0, 0.48, 0.47, 426.0, 40000]]) y_pred = model.predict(x_new)

Code 4: Code for training the Model Predictive Control (MPC) model with KKT-HardNet.

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. IPOPT PINN DC3 KKT-HardNet

Objective

Mean eq.

Max eq.

Mean ineq.

Max ineq.

1.32 × 10−1

4.75 × 10−10

5.63 × 10−8

4.03 × 10−14

9.98 × 10−9

1.19 × 10−1

5.31 × 10−3

8.16 × 10−2

6.18 × 10−6

2.72 × 10−2

Training Time (s)

Inference Time (ms/batch)

−− 137.68 2.70 × 10−1 0.00 0.00 1.39 × 103 3.45 × 104 1821.28 −1 −5 −3 1.60 × 10 2.28 × 10 1.48 × 10 0.00 0.00 10364.44 *Solver time is reported as average solve time considering full parallelization across 100 instances

6.30∗ 0.72 31.60 166.07

As shown in Table 8, KKT-HardNet obtains an average objective value of 1.60×10−1 compared to 1.32×10−1 obtained by the reference optimizer, while maintaining zero inequality violation and small equality residuals. The PINN obtains a lower objective value but remains infeasible, whereas DC3 satisfies the equality constraints but exhibits large inequality violations. The reported inference times correspond to batches of 100 instances; however, closed-loop control requires sequential single-instance evaluations, which are analyzed separately below. We integrate the previous RTO model, which finds the optimal set points, with the MPC model to complete the closed loop control, as illustrated in Figure 12. The economically optimal operating targets obtained from the RTO layer are used as set points for the MPC. The closed-loop RTO–MPC response is shown in Fig. 13. The RTO layer is updated every 1000 s, while the MPC is evaluated every 5 s. A perfect actuator is assumed [39], and the nonlinear CSTR is propagated using fourth-order Runge–Kutta integration. The updated state is supplied to the controller at the next sampling instance, forming the receding-horizon closed loop. Figure 13 and Table 9 illustrate that the controller closely tracks the economically optimal set points. The RMSE values for 𝐶 𝐴, 𝐶 𝐵 , and 𝑇 are 6.35 × 10−4 mol/L, 6.34 × 10−4 mol/L, and 0.278 K, respectively. Following each RTO update, the process settles within approximately 185–255 s, well within the 1000 s RTO interval, with maximum overshoot below 0.03% for all controlled variables. Constraint satisfaction was verified over all 2000 controller evaluations. The mean and maximum absolute equality residuals were 1.90 × 10−9 and 5.15 × 10−7 , respectively, while the mean and maximum inequality violations were 7.72 × 10−13 and 2.95 × 10−9 . No violations were observed for the predicted state bounds, heat-input bounds, horizon control-move constraints, or simulated plant operating bounds. The remaining small inequality residuals arise from the Lyapunov condition and remain below the 10−8 projection tolerance. 24

Heat Price (p)

KKTHN-RTO

x sp , Qsp

KKTHN-MPC

Q

CSTR System x

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 (𝐶 𝐴𝑠 𝑝 , 𝐶 𝐵𝑠 𝑝 , 𝑇 𝑠 𝑝 , 𝑄 𝑠 𝑝 ) , while the KKT-HardNet MPC regulates the nonlinear CSTR toward the updated operating conditions.

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. Variable 𝐶𝐴 𝐶𝐵 𝑇

RMSE

Mean settling time (s)

Max. overshoot (%)

Mean steady-state abs. error

6.35 × 10−4 6.34 × 10−4 0.278

222.8 223.3 204.4

0.026 0.014 0.004

5.27 × 10−7 2.85 × 10−7 3.47 × 10−5

25

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.

The Lyapunov based analysis in Fig. 14 shows that 𝑉 ( 𝑥) ¯ increases immediately after each reference update and subsequently decays toward zero. Across all 2000 evaluations, the condition 𝑉¤ ( 𝑥, ¯ 𝑢 0 ) ≤ 𝑉¤ ( 𝑥, ¯ ℎ( 𝑥)) ¯ −8 −9 exhibits no violations above the 10 tolerance. The minimum Lyapunov margin is −2.95 × 10 and the mean margin is 1.19 × 10−5 , indicating consistent satisfaction of the imposed Lyapunov condition. This behavior shows that the controller drives the system toward the updated equilibrium after each set-point change, providing numerical evidence of closed-loop convergence and stability. The consistent satisfaction of the Lyapunov condition further supports the stabilizing behavior of the controller. The real-time computational performance of the MPC was evaluated over 2000 warm single-instance predictions. The mean and median inference latencies were 20.20 and 19.94 ms, while the 95th and 99th percentile latencies were 29.92 and 34.86 ms. The maximum latency was 79.41 ms, corresponding to only 1.59% of the 5 s sampling interval, and no deadline misses were observed. The first evaluation required approximately 35.95 s because of one-time JAX compilation and was excluded from the online timing statistics. At the nine instances where both RTO and MPC were evaluated sequentially, the combined median and maximum latencies were 50.95 and 59.48 ms, respectively. The closed-loop economic performance is evaluated using 𝐽plant = 𝐶 𝐴/𝐶 𝐴0 +𝑝𝑄 and 𝐽RTO = 𝐶 𝐴𝑠 𝑝 /𝐶 𝐴0 +𝑝𝑄 𝑠 𝑝 . Over the complete 10,000 s simulation, the mean plant and RTO objectives were 0.54018 and 0.54025, respectively, with a mean absolute economic gap of 2.12 × 10−4 . The maximum absolute gap was 4.50 × 10−3 and occurred immediately after an RTO update, after which the MPC rapidly returned the process objective toward the RTO target. Additional analyses of control effort, robustness, alternative heat-price profiles, and local controller sensitivity are provided in Appendix D. 26

Figure 15: Closed-loop economic performance of the coupled RTO–MPC framework. The top panel compares the instantaneous plant objective 𝐽plant with the steady-state RTO target 𝐽RTO , while the remaining panels show the instantaneous and cumulative economic gaps.

6

PCML for Grey-Box Optimization

Data-driven optimization is useful when the relationship between operating conditions and process performance is available primarily through experiments, sensor measurements, or high-fidelity simulations rather than an explicit differentiable process model. A surrogate model can first be trained to approximate this input–output relationship and subsequently used as an efficient model for optimization. In general, let 𝒙 ∈ X denote the controllable operating conditions and 𝒚 = M(𝒙) the outputs predicted by a trained surrogate M. The downstream optimization problem can then be expressed as 𝒙 ∗ = arg min 𝒙∈X 𝐹 (𝒙, M(𝒙)) , where 𝐹 represents the desired process-performance or economic objective. The reliability of this approach depends on the accuracy and physical consistency of the learned surrogate, since the optimizer may otherwise identify operating conditions that exploit prediction errors or produce physically inconsistent outputs. KKT-HardNet provides a differentiable surrogate whose predictions satisfy the physical constraints incorporated during training. Consequently, gradients of the downstream objective can be propagated through the trained KKT-HardNet model to optimize the input conditions directly, while the surrogate predictions remain constrained by the available domain knowledge. This provides a general framework for using observation-based models not only for prediction, but also for identifying improved operating conditions. 6.1

Simulation-based Optimization of a Syngas Production Process

We consider the syngas production pathway illustrated in Fig. 16. The process uses CO2 and H2 as feed streams, both entering the mixer at a temperature of 25◦ C and a pressure of 1 atm. The CO2 feed flow rate is fixed at 1 kmol, h−1 as the basis, while the H2 feed flow rate is allowed to vary. We define the ratio between the CO2 and H2 feed flow rates as 𝑎 = 𝐹CO2 /𝐹H2 , where 𝐹CO2 and 𝐹H2 denote the molar flow rates of CO2 and H2 , respectively. 27

COMP-101

CO2

Syngas

HE-101

VAL-101 REAC-101

HE-102

H2

SEP-101 MIX-101

Water

Figure 16: A syngas production pathway The two feed streams are first combined in a mixer, after which the mixed gas stream is compressed from the feed pressure to the desired operating pressure, 𝑃op , using an isentropic compressor. The compressor is driven by an electrical utility. Following compression, a simple heat exchanger is used to bring the mixed gas stream to the desired reactor operating temperature, 𝑇op . Depending on the temperature of the compressed stream relative to the specified reactor operating temperature, the heat exchanger can provide either heating or cooling. The heating requirement can be supplied using a heating utility, such as steam or electricity, whereas the cooling requirement can be supplied using a cooling utility, such as cooling water or refrigerant. In this manner, the feed stream entering the reactor is maintained at the specified operating temperature 𝑇op and operating pressure 𝑃op . The conditioned mixed-gas stream is subsequently introduced into a Gibbs equilibrium reactor. Instead of employing an explicit kinetic reaction model, the equilibrium reactor determines the product composition by minimizing the total Gibbs free energy of the system subject to elemental conservation constraints. The chemistry occurring in the reactor can primarily be described by the CO2 methanation and reverse water–gas shift reactions, CO2 + 4H2 ⇌ CH4 + 2H2 O,

(16)

CO2 + H2 ⇌ CO + H2 O.

(17)

The Gibbs equilibrium reactor is configured with a maximum of 50 equilibrium iterations. The reactor therefore predicts the equilibrium distribution among the species CO2 , H2 , CO, CH4 , and H2 O at the specified operating conditions 𝑇op and 𝑃op . Following the reactor, the product stream is passed through a valve, where it is depressurized from the reactor operating pressure to approximately 1 atm. The depressurized stream is subsequently passed through a cooler, where its temperature is reduced to 50◦ C. The cooled reactor effluent is then introduced into a vapor–liquid separator operating at approximately 𝑇sep = 25◦ C and 1 atm. The separator performs vapor–liquid equilibrium separation and produces two outlet streams. The top stream corresponds to the gas-phase product, referred to as the syngas stream, and contains five gas-phase components: CO2 , H2 , CO, CH4 , and H2 O. The bottom stream corresponds to the liquid-phase product and primarily contains condensed water. The Aspen Plus simulation employs the Peng–Robinson equation of state as the thermodynamic property method for evaluating the required phase and equilibrium properties throughout the flowsheet. To generate the dataset used in this work, we sample 4000 operating points by varying the feed-flow ratio 𝑎 , reactor operating temperature 𝑇op , and reactor operating pressure 𝑃op . The operating points are generated using Latin hypercube sampling with a random seed of 42. The corresponding sampling ranges are 𝑎 ∈ [1, 10] , 𝑇op ∈ [100, 1000] ◦ C, and 𝑃op ∈ [1, 20] atm. For each sampled combination of 𝑎 , 𝑇op , and 𝑃op , the complete Aspen Plus flowsheet is simulated to obtain the corresponding equilibrium product composition and process utility requirements. Prediction Model: The syngas formulation has 3 parameters, 6 primal variables, 3 equality constraints, and no inequality constraints. The input parameters are the feed ratio 𝑎 , operating pressure 𝑃op , and operating 28

temperature 𝑇op , while the predicted variables are the outlet molar flows of CO2 , H2 , CO, CH4 , and H2 O, together with the total utility cost. The three equality constraints represent the elemental carbon, oxygen, and hydrogen balances and are given as 𝐹CO2 + 𝐹CO + 𝐹CH4 = 1 2𝐹CO2 + 𝐹CO + 𝐹H2 O = 2 (18) 4𝐹CH4 + 2𝐹H2 + 2𝐹H2 O = 2𝑎. We use the 4000 samples generated from the Aspen Plus simulations and split the dataset in an 80 : 20 ratio for training and validation respectively. The network structure consists of 2 hidden layers with 64 hidden nodes in each layer and ReLU activation. Since the formulation contains only linear equality constraints, we consider the “linear-equality" projection formulation [15], where the network predicts the 6 primal variables together with 3 equality dual variables. We first perform MLP pretraining for 100 epochs and subsequently activate the hard-constraint projection. The model is trained for a total of 2000 epochs with a batch size of 50, learning rate of 10−3 , random seed 42. The training loss consists of the mean-squared data loss together with consistency losses for the primal and equality-dual variables, for which both consistency weights are set to 1. Since no inequality constraints are present in this formulation, no inequality-dual variables are required. The linear-equality projection enforces the three elemental balance constraints during the projection stage, resulting in equality-constraint residuals close to machine precision. The code to train the model is given in codeblock 5. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34

import os os.environ["JAX_ENABLE_X64"] = "True" from pathlib import Path import numpy as np import jax jax.config.update("jax_enable_x64", True) import jax.numpy as jnp from kkthn import KKTHardNet, Loss, MLP, Projection SEED = 42 TRAIN = { "epochs": 2000,"epoch_mlp": 100,"batch_size": 50,"learning_rate": 1e-3,"train_frac": 0.8," seed": SEED,"dtype": "float64","print_every": 100,"minimum_projection_convergence": 0.95," readiness_extra_epochs": 500,"readiness_check_every": 10,"backend": "cpu","output_dir": str( Path.cwd() / "runs"),} DATA_PATH = Path.cwd() / "syngas" / "syngas_dataset_lhc.csv" # Backbone nn = MLP(hidden_sizes=(64, 64), activation="relu", predict_duals=True, inequality_dual_floor=1e -12) # Projection projection = Projection(name="syngas_production", method="linear-equality") x = projection.add_parameter(["a", "P", "T"]) y = projection.add_variable(["CO2", "H2", "CO", "CH4", "H2O", "Util"]) projection.constraints.add( y.CO2 + y.CO + y.CH4 == 1, 2.0 * y.CO2 + y.CO + y.H2O == 2, 4.0 * y.CH4 + 2.0 * y.H2 + 2.0 * y.H2O == 2.0 * x.a ) # Loss loss = Loss.mse(primal_consistency_weight=1.0, equality_dual_consistency_weight=1.0, inequality_dual_consistency_weight=1.0)

29

35 36 37 38 39 40 41 42 43

# Model Definition model = KKTHardNet("surrogate", nn, loss, projection, train=TRAIN) model.dataset(parameters=DATA_PATH, variables=DATA_PATH) # Training and Summary linear_equality_result = model.model() model.summary() model.plot_history()

Code 5: Code for training the syngas model with KKT-HardNet.

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. MLP PINN ENFORCE KKT-HardNet

MSE

Mean eq.

Max eq.

Mean ineq.

Max ineq.

Training Time (s)

Inference Time (ms/batch)

6.7 × 10−5 1.31 × 10−4 3.45 × 10−5 4.48 × 10−5

5.40 × 10−3 6.47 × 10−3 1.31 × 10−7 2.25 × 10−16

0.11 0.05 3.81 × 10−6 3.55 × 10−15

– – – –

– – – –

92.02 90.87 226.70 60.80

0.2 0.2 0.66 1.73

Table 10 compares the predictive accuracy, elemental-balance satisfaction, training cost, and inference performance of the surrogate models for the syngas production case study. KKT-HardNet achieves an MSE of 4.48 × 10−5 while satisfying the three elemental-balance equality constraints to essentially machine precision, with mean and maximum residuals of 2.25 × 10−16 and 3.55 × 10−15 , respectively. The unconstrained MLP and PINN attain similar prediction errors but exhibit substantially larger equality residuals, whereas ENFORCE also enforces the balances accurately but requires a longer training time. KKT-HardNet therefore combines high predictive accuracy with exact elemental-balance consistency and the lowest training time among the hard-constrained approaches considered, although its inference time is higher than those of the unconstrained and soft-constrained neural-network baselines. Optimal operating conditions: After training, the KKT-HardNet model is used as a differentiable surrogate to optimize the syngas operating conditions with respect to the feed ratio 𝑎 , operating pressure 𝑃op , and operating temperature 𝑇op . The scalarized objective is defined as 𝐹 = 𝛼𝑈 − (1 − 𝛼) 𝑋CO2 , where 𝑋CO2 = 1 − 𝐹CO2 ,out and 𝛼 ∈ [0, 1] controls the trade-off between utility consumption and CO2 conversion. The optimization is performed using a multi-start projected Adam scheme in the normalized input space, as shown in Figure 17, with the operating variables constrained to their prescribed bounds. To characterize the trade-off in greater detail, the problem is solved for 𝛼 = 0, 0.02, . . . , 1.0 using 100 random initializations for each value of 𝛼. The resulting operating conditions and performance metrics are shown in Fig. 18. For small 𝛼, the optimizer emphasizes conversion, resulting in operating points with nearly complete CO2 conversion and relatively high utility consumption. Over the broad interval 𝛼 ≈ 0.06–0.74, the optimizer repeatedly identifies an operating condition near 𝑎 ≈ 4.79, 𝑃op = 1 atm, and 𝑇op = 100 ◦ C, for which the predicted conversion remains close to 0.999 and the utility is approximately 0.345. A pronounced transition occurs near 𝛼 = 0.76, after which the solution shifts toward lower feed ratio and lower utility consumption, accompanied by a substantial reduction in CO2 conversion. For 𝛼 ≳ 0.86, the optimum reaches the lower operating bounds 𝑎 = 1, 𝑃op = 1 atm, and 𝑇op = 100 ◦ C, with a predicted utility of approximately 0.110 and conversion of approximately 0.257. Since the KKT-HardNet surrogate enforces the elemental-balance constraints by construction, the outer optimization is carried out directly over the bounded operating variables without explicitly reintroducing these balance equations. 30

Trained surrogate M from KKT-HardNet and optimization parameters Select 𝛼 and set 𝐹 ∗ ← +∞ Sample (𝑎, 𝑃op , 𝑇op ) and initialize Adam Evaluate M and objective 𝐹 = 𝛼𝑈 − (1 − 𝛼) 𝑋CO2 Compute gradient, Adam step, and projection Next iteration

Converged or 𝑘 = 𝐾?

No

Next start

Next 𝛼

Yes

If 𝐹 < 𝐹 ∗ , update and store best solution

More random starts?

Yes

No

Yes

More 𝛼 values? No

Return optimal operating conditions for each 𝛼

Figure 17: Multi-start projected Adam for syngas operating-condition optimization. More details about the algorithm is given in Appendix E

31

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 𝑎 , operating pressure 𝑃op , and operating temperature 𝑇op .

7

Conclusions

Integrating optimization algorithms with machine learning provides a promising route for enforcing physical and operational constraints in data-driven models. In this work, we developed an efficient, scalable, and open-source implementation of KKT-HardNet with modular support for custom neural-network backbones and projection algorithms. Implicit differentiation enables efficient training through iterative projection layers without unrolling the solver iterations. Across applications in process simulation, real-time optimization and control, and data-driven optimization, KKT-HardNet consistently achieved strong predictive or optimization performance while maintaining tight constraint satisfaction. The results also demonstrate the flexibility of the framework across different problem structures and projection strategies. Future work will focus on larger-scale applications and problem-specific projection methods that further exploit sparsity and structure. More broadly, we aim to establish KKT-HardNet as a practical development and benchmarking platform for physics-constrained machine learning.

Acknowledgment The authors gratefully acknowledge partial funding support from Texas A&M University.

Declaration of Generative AI and AI-Assisted Technologies The authors used generative artificial intelligence (AI) tools to assist with generating initial code, language editing, grammar, readability, and clarity. Generative AI tools were not used to generate research data, perform scientific analyses, develop conclusions, or make scientific claims. All AI-assisted contents were reviewed, tested and modified by the authors, who take full responsibility for the scientific content and integrity of the manuscript. 32

References [1] L. Von Rueden, S. Mayer, K. Beckh, B. Georgiev, S. Giesselbach, R. Heese, B. Kirsch, J. Pfrommer, A. Pick, R. Ramamurthy, et al., “Informed machine learning–a taxonomy and survey of integrating prior knowledge into learning systems,” IEEE Transactions on Knowledge and Data Engineering, vol. 35, no. 1, pp. 614–633, 2021. [2] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang, “Physics-informed machine learning,” Nature Reviews Physics, vol. 3, no. 6, pp. 422–440, 2021. [3] M. Raissi, P. Perdikaris, and G. E. Karniadakis, “Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations,” Journal of Computational physics, vol. 378, pp. 686–707, 2019. [4] S. Cai, Z. Wang, S. Wang, P. Perdikaris, and G. E. Karniadakis, “Physics-informed neural networks for heat transfer problems,” Journal of Heat Transfer, vol. 143, no. 6, p. 060801, 2021. [5] Z. Mao, A. D. Jagtap, and G. E. Karniadakis, “Physics-informed neural networks for high-speed flows,” Computer Methods in Applied Mechanics and Engineering, vol. 360, p. 112789, 2020. [6] D. W. Abueidda, P. Pantidis, and M. E. Mobasher, “Deepokan: Deep operator network based on kolmogorov arnold networks for mechanics problems,” Computer Methods in Applied Mechanics and Engineering, vol. 436, p. 117699, 2025. [7] A. Yazdani, L. Lu, M. Raissi, and G. E. Karniadakis, “Systems biology informed deep learning for inferring parameters and hidden dynamics,” PLoS computational biology, vol. 16, no. 11, p. e1007575, 2020. [8] N. Ahmadi Daryakenari, M. De Florio, K. Shukla, and G. E. Karniadakis, “AI-Aristotle: A physicsinformed framework for systems biology gray-box identification,” PLOS Computational Biology, vol. 20, no. 3, p. e1011916, 2024. [9] Y. Zhu, N. Zabaras, P.-S. Koutsourelakis, and P. Perdikaris, “Physics-constrained deep learning for high-dimensional surrogate modeling and uncertainty quantification without labeled data,” Journal of computational physics, vol. 394, pp. 56–81, 2019. [10] A. Iftakher, R. Golder, B. N. Roy, and M. M. F. Hasan, “Physics-informed neural networks with hard nonlinear equality and inequality constraints,” Computers & Chemical Engineering, p. 109418, 2025. [11] A. Mukherjee and V. M. Zavala, “Physics-constrained machine learning for chemical engineering,” Current Opinion in Chemical Engineering, vol. 51, p. 101228, 2026. [12] T. Beucler, M. Pritchard, S. Rasp, J. Ott, P. Baldi, and P. Gentine, “Enforcing analytic constraints in neural networks emulating physical systems,” Physical review letters, vol. 126, no. 9, p. 098302, 2021. [13] A. S. Zamzam and K. Baker, “Learning optimal solutions for extremely fast ac optimal power flow,” in 2020 IEEE international conference on communications, control, and computing technologies for smart grids (SmartGridComm), pp. 1–6, IEEE, 2020. [14] P. Donti, D. Rolnick, and J. Z. Kolter, “DC3: A learning method for optimization with hard constraints,” in International Conference on Learning Representations, 2021. [15] H. Chen, G. E. C. Flores, and C. Li, “Physics-informed neural networks with hard linear equality constraints,” Computers & Chemical Engineering, vol. 189, p. 108764, 2024. [16] G. E. C. Flores, H. Chen, and C. Li, “Enforcing hard linear constraints in deep learning models with decision rules,” in The Thirty-ninth Annual Conference on Neural Information Processing Systems, 2025. [17] Y. Min and N. Azizan, “Hardnet: Hard-constrained neural networks with universal approximation guarantees,” arXiv preprint arXiv:2410.10807, 2024. 33

[18] P. D. Grontas, A. Terpin, E. C. Balta, R. D’Andrea, and J. Lygeros, “Pinet: Optimizing hard-constrained neural networks with orthogonal projection layers,” arXiv preprint arXiv:2508.10480, 2025. [19] G. Lastrucci and A. M. Schweidtmann, “ENFORCE: Nonlinear constrained learning with adaptivedepth neural projection,” arXiv preprint arXiv:2502.06774, 2025. [20] B. N. Roy, R. Golder, and M. Hasan, “NLPOpt-Net: A learning method for nonlinear optimization with feasibility guarantees,” arXiv preprint arXiv:2605.00260, 2026. [21] A. Goertzen, K. Alim, Y. Min, and N. Azizan, “Hardnet++: Nonlinear constraint enforcement in neural networks,” IEEE Control Systems Letters, 2026. [22] F. M. Mohammadi, H. Budman, and J. L. Pulsipher, “PL-KKT-hPINN: Enforcing nonlinear equality constraints on neural networks via piecewise-linear projection,” arXiv preprint arXiv:2606.10682, 2026. [23] E. N. Pistikopoulos, N. A. Diangelakis, and R. Oberdieck, Multiparametric optimization and control. John Wiley & Sons, 2020. [24] K. Gurney, An introduction to neural networks. CRC press, 2018. [25] P. Shah, S. Pahari, R. Bhavsar, and J. S.-I. Kwon, “Hybrid modeling of first-principles and machine learning: A step-by-step tutorial review for practical implementation,” Computers & Chemical Engineering, vol. 194, p. 108926, 2025. [26] G. Cybenko, “Approximation by superpositions of a sigmoidal function,” Mathematics of control, signals and systems, vol. 2, no. 4, pp. 303–314, 1989. [27] S. Haykin, Neural networks and learning machines, 3/E. Pearson education india, 2009. [28] A. D. Rasamoelina, F. Adjailia, and P. Sinčák, “A review of activation function for artificial neural network,” in 2020 IEEE 18th world symposium on applied machine intelligence and informatics (SAMI), pp. 281–286, IEEE, 2020. [29] J. Nocedal and S. J. Wright, Numerical optimization. Springer, 2006. [30] B. Amos and J. Z. Kolter, “OptNet: Differentiable optimization as a layer in neural networks,” in International conference on machine learning, pp. 136–145, PMLR, 2017. [31] R. Golder, B. N. Roy, and M. M. F. Hasan, “Physics constrained machine learning for modeling and optimization of chemical process systems,” Systems and Control Transactions, vol. 5, pp. 1925–1931, 2026. [32] D. E. Rumelhart, G. E. Hinton, and R. J. Williams, “Learning representations by back-propagating errors,” nature, vol. 323, no. 6088, pp. 533–536, 1986. [33] M. R. Osborne, “Fisher’s method of scoring,” International Statistical Review/Revue Internationale de Statistique, pp. 99–117, 1992. [34] K. Levenberg, “A method for the solution of certain non-linear problems in least squares,” Quarterly of applied mathematics, vol. 2, no. 2, pp. 164–168, 1944. [35] D. W. Marquardt, “An algorithm for least-squares estimation of nonlinear parameters,” Journal of the society for Industrial and Applied Mathematics, vol. 11, no. 2, pp. 431–441, 1963. [36] J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, Y. Katariya, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, and Q. Zhang, “JAX: composable transformations of Python+NumPy programs,” 2018. [37] A. Wächter and L. T. Biegler, “On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming,” Mathematical programming, vol. 106, no. 1, pp. 25–57, 2006. [38] A. Kumar and P. Daoutidis, Control of nonlinear differential algebraic equation systems with applications to chemical processes, vol. 397. CRC Press, 1999. 34

[39] Z. Zhang, Z. Wu, D. Rincon, and P. D. Christofides, “Real-time optimization and control of nonlinear processes using machine learning,” Mathematics, vol. 7, no. 10, p. 890, 2019.

35

A

Benchmarking: metrics and other methods

We compare the performance of KKT-HardNet with other existing methods on the following metrics. Data Loss: Data loss refers to how close the final prediction is to the available data. We measure the data loss using the standard Mean Squared Error (MSE) metric. Optimality: Optimality refers to how close the final prediction is to the true optimal solution. This is applicable to all the case studies on optimization. Feasibility: Feasibility is measured by the constraint violation. The constraint violation is computed as the absolute residual i.e, |ℎ(𝑥, 𝑦)| for equality constraint, and 𝑔+ (𝑥, 𝑦) for inequality where 𝑔+ = max(0, 𝑔) . For a feasible point, both the equality and inequality violations should be zero. We evaluate both the maximum violation and the mean violation. Speed: Speed refers to the time required for inference after the training is finished. The inference time is given for each batch. For optimizer we consider full parallelization and for a fare comparison we divide the total solve time by total number of instances in a batch. For inference, batch size across the case studies remains the same which is 100 instances. We compare KKT-HardNet against the following methods: • Optimizer: A traditional optimization solver. We consider IPOPT [37] as the state-of-the-art optimizer. • NN: A regular artificial neural network model with no soft loss considered. We compare with MLP for the surrogate modeling case studies. • Soft Constrained NN: We use Physics Informed Neural Networks (PINNs) [3] as the soft constrained methods to compare our model with. • ENFORCE: ENFORCE [19] employs an adaptive depth iterative projection method on the approximated (linearized) constraints. We use the github codespace for ENFORCE including default tunable parameters for comparison. • DC3: The deep constraint completion and correction (DC3) [14] method is state-of-the-art that uses predict and completion method to satisfy the equality constraints. The inequality correction is then done by mapping infeasible points to feasible points using an internal gradient descent step during training. We use the github codespace for DC3 including default tunable parameters for comparison. Across all the case studies we present, a significant portion of the code remains same. These include library imports, backbone setup, training configurations, and solver configurations.

B

Methods Included to Solve the Projection Problem

We briefly describe reformulation for the KKT conditions to find a feasible solution, a general algorithm to solve the projection optimization problem and the explicit projection for linear equality that we include in our implementation. B.1

Explicit Projection: Linear Equality Constraints

When the constraint space is defined with linear equalities of the form 𝐴𝒙 + 𝐵𝒚 = 𝑏 the KKT condition of the projection in Eq. 4 yields the analytical solution: 𝒚˜ = 𝒚ˆ − 𝐵⊤ (𝐵𝐵⊤ ) −1 (𝐵 𝒚ˆ + 𝐴𝒙 − 𝑏),

(19)

which can be written in affine form 𝒚˜ = 𝐴∗ 𝒙 + 𝐵∗ 𝒚ˆ + 𝑏 ∗ where 𝐴∗ = 𝐵⊤ (𝐵𝐵⊤ ) −1 𝐴, 𝐵∗ = 𝐼 − 𝐵⊤ (𝐵𝐵⊤ ) −1 𝐵, and 𝑏 ∗ = 𝐵⊤ (𝐵𝐵⊤ ) −1 𝑏 (see [15] for details). This explicit affine form of the solution allows to have a 36

projection with single iteration to ensure feasibility. 𝐴∗ , 𝐵∗ and 𝑏 ∗ can be computed once and used during training and inference making the projection computationally inexpensive. B.2

Implicit Projection I: KKT conditions via Fischer-Burmeister Reformulation

With 𝝀 ∈ R 𝑝 and 𝝁 ∈ R𝑞 being the Lagrangian multipliers for the equality and inequality constraints, the KKT conditions of the projection in Eq. 4 for general nonlinear constraints can be written as the following square system of nonlinear equation with slack variables 𝒔 ∈ R𝑞 corresponding to the inequality constraints (see [10] for details). 𝒚 − 𝒚ˆ +

𝑝 ∑︁

𝜆𝑖 ∇𝒚 ℎ𝑖 (𝒙, 𝒚) +

𝑖=1

𝑞 ∑︁

𝜇 𝑗 ∇𝒚 𝑔 𝑗 (𝒙, 𝒚) = 0,

(20a)

𝑗=1

ℎ𝑖 (𝒙, 𝒚) = 0, 𝑔 𝑗 (𝒙, 𝒚) + 𝑠 𝑗 = 0, √︃ 𝜇 𝑗 + 𝑠 𝑗 − 𝜇2𝑗 + 𝑠2𝑗 = 0,

𝑖 = 1, 2, . . . , 𝑝 𝑗 = 1, 2, . . . , 𝑞

(20b) (20c)

𝑗 = 1, 2, . . . , 𝑞.,

(20d)

where, Eq. 20a is the stationarity, 20b and 20c are primal feasibility conditions, and Eq. 20d is the FischerBurmeister reformulation of the dual and complementarity conditions without losing any fidelity. Eq. 20 can be solved with an iterative Gauss-Newton procedure with Armijo backtracking to find a stationary solution i.e. feasible by construction. It is trivial to show that, if the constraints are only linear equality then Eq. 20 yields the same expression as Eq. 19. B.3

Implicit Projection II: Primal-Dual Interior Point Method

Check [37] for details on the theory and algorithms. We adopt a primitive implementation of the algorithm to fit for batch inference and suitable for GPU runs.

C

Reduction and Heuristics on KKT Conditions via Fishcer-Burmeister Reformulation

In order to solve the KKT system given in Equation 20, a damped Newton method can be applied. We consider 𝑦 𝜇 (𝑘) as the iteration counter given an iterate ( 𝑦 (𝑘 ) , 𝜆 (𝑘 ) , 𝜇 (𝑘 ) , 𝑠 (𝑘 ) ). The search direction (𝑑 (𝑘 ) , 𝑑 𝜆(𝑘 ) , 𝑑 (𝑘 ) , 𝑑 𝑠(𝑘 ) ) are obtained from the following linearization: 𝐻 (𝑘 )  ⊤  𝐴 (𝑘 )  ⊤ 𝐶  (𝑘 )  0 

𝐴 (𝑘 ) 0 0 0

𝐶 (𝑘 ) 0 0 𝜇 𝐺 (𝑘 )

𝑦  𝑑 (𝑘 )   𝑦 (𝑘 ) − 𝑦ˆ + 𝐴 (𝑘 ) 𝜆 (𝑘 ) + 𝐶 (𝑘 ) 𝜇 (𝑘 )   𝜆     𝑑 (𝑘 )    ℎ(𝑥, 𝑦 (𝑘 ) )   𝜇  =− .  𝑑    𝑔(𝑥, 𝑦 (𝑘 ) ) + 𝑠 (𝑘 )   (𝑘 )    𝑠    𝑠 Φ(𝜇 , 𝑠 ) 𝐺 (𝑘 )  𝑑 (𝑘 )  (𝑘 ) (𝑘 )    

0 0 𝐼

(21)

Here, Φ(𝜇, 𝑠) is the Fischer-Burmeister equation block, 𝐻 (𝑘 ) := 𝐼 + ∇2𝑦 𝑦 [𝜆⊤ ℎ(𝑥, 𝑦) + 𝜇⊤ 𝑔(𝑥, 𝑦)] is the 𝜇 Hessian of the Lagrangian, 𝐴 (𝑘 ) := ∇ 𝑦 ℎ(𝑥, 𝑦) , 𝐶 (𝑘 ) := ∇ 𝑦 𝑔(𝑥, 𝑦) , 𝐺 (𝑘 ) = ∇ 𝜇 Φ and 𝐺 𝑠(𝑘 ) = ∇𝑠 Φ. Note that, 𝜇 both 𝐺 (𝑘 ) and 𝐺 𝑠(𝑘 ) are diagonal matrix defined as follows: 𝜇

𝐺 (𝑘 ) = diag (𝜙 𝜇 𝑗 ); 𝐺 𝑠(𝑘 ) = diag (𝜙 𝑠 𝑗 ); 𝜙 𝜇 𝑗 = 1 − √︃

𝜇𝑗 𝜇2𝑗 + 𝑠2𝑗

; 𝜙 𝑠 𝑗 = 1 − √︃

𝑠𝑗 .

(22)

𝜇2𝑗 + 𝑠2𝑗

While the gradients in Equation 22 are not defined for (0, 0) , a smooth Fischer-Burmeister formulation can be used to get the gradient in such cases. Instead of solving the nonsymmetric linear system shown in Equation 21 directly, following smaller symmetric system can be solved first, 37

𝐻 (𝑘 )  ⊤  𝐴 (𝑘 )  ⊤ 𝐶  (𝑘 )

𝑦    𝑑 (𝑘 )  𝐶 (𝑘 ) 𝑦 (𝑘 ) − 𝑦ˆ + 𝐴 (𝑘 ) 𝜆 (𝑘 ) + 𝐶 (𝑘 ) 𝜇 (𝑘 )    𝜆  ,  𝑑  = −  0 ℎ(𝑥, 𝑦 (𝑘 ) ) (𝑘 )    −1 𝑠 −1 𝜇   𝜇  𝑠  𝑔(𝑥, 𝑦 (𝑘 ) ) + 𝑠 (𝑘 ) − 𝐺 (𝑘 ) Φ(𝜇 (𝑘 ) , 𝑠 (𝑘 ) )  −𝐺 (𝑘 ) 𝐺 (𝑘 )  𝑑 (𝑘 )    

𝐴 (𝑘 ) 0 0

(23)

with eliminating the last block. The vector 𝑑 𝑠(𝑘 ) can be then obtained from 𝑑 𝑠(𝑘 ) = 𝐺 𝑠(𝑘 ) −1 Φ(𝜇 (𝑘 ) , 𝑠 (𝑘 ) ) − 𝐺 𝑠(𝑘 ) −1 𝐺 (𝑘 ) 𝑑 (𝑘 ) . 𝜇

𝜇

(24)

While, to hold this true 𝐺 𝑠(𝑘 ) has to be invertible, which is often time not at the solution, we bypass this with some elementwise diagonal pivoting into the 𝐺 𝑠(𝑘 ) matrix retaining the reduced linear system symmetric. Furthermore, if 𝐴 (𝑘 ) is not full rank the system will be singular and a solution might not exist, therefore regularization is used while solving this matrix. After the search direction is obtained the next iterate is obtained with 𝑦 𝑦 𝑦 (𝑘+1) = 𝑦 (𝑘 ) + 𝛼 𝑘 𝑑 (𝑘 ) 𝜆 (𝑘+1) = 𝜆 (𝑘 ) + 𝛼𝜆𝑘 𝑑 𝜆(𝑘 ) 𝜇 𝜇

𝜇 (𝑘+1) = 𝜇 (𝑘 ) + 𝛼 𝑘 𝑑 (𝑘 )

(25)

𝑠 (𝑘+1) = 𝑠 (𝑘 ) + 𝛼 𝑘𝑠 𝑑 𝑠(𝑘 ) 𝑦

𝜇

for some 𝛼 (𝑘 ) , 𝛼𝜆(𝑘 ) , 𝛼 (𝑘 ) , 𝛼 𝑠(𝑘 ) ∈ (0, 1] . In our implementation we kept both methods ‘fischer-burmeister-newton’ and ‘fischer-burmeister-reduced’ with the full linear system and reduced linear system respectively. Moreover, the reduction is only on the slack variables. Hence the performance change is only visible when the constraint set includes inequalities.

D

Additional Closed-Loop Analysis of the RTO–LMPC Framework

D.1

LMPC Control-Effort Analysis

The LMPC control action is expressed in deviation form as 𝑢 0 = 𝑄 − 𝑄 𝑠 𝑝 , allowing the transient corrective action of the controller to be separated from the economically optimal steady-state heat input supplied by the RTO layer. Figure 19 shows that the applied heat input remains very close to 𝑄 𝑠 𝑝 throughout the closed-loop simulation, with appreciable corrective action occurring primarily immediately after an RTO set-point update. The mean absolute control correction was only 0.736 cal/s, with an RMS value of 2.14 cal/s and a maximum magnitude of 13.81 cal/s. This is small relative to the nominal heating rate, which varied between approximately 2.34 × 104 and 3.96 × 104 cal/s. The total variation of the applied heat input was 2.34 × 104 cal/s, compared with 2.32 × 104 cal/s for the RTO steady-state heat-input trajectory. Thus, the LMPC correction increased the total control variation by only 198.8 cal/s, or approximately 0.86%, indicating that most of the overall input variation originates from changes in the economic operating target rather than aggressive feedback action. The largest observed change between successive implemented heat inputs was 3.26 × 103 cal/s, and the applied heat input remained well within its physical bounds throughout the simulation. D.2

Closed-Loop Performance under Challenging Heat-Price Profiles

To evaluate the closed-loop response beyond the nominal staircase heat-price trajectory, the coupled RTO– LMPC framework was tested using sinusoidal, random piecewise-constant, sharp-spike, and monotonic-ramp price profiles. Figure 20 shows that the controller remains stable and bounded for all investigated cases, with zero physical plant-bound violation in every simulation. The baseline and ramp profiles resulted in the smallest tracking errors, with temperature RMSE values of 0.278 K and 0.299 K, respectively. More 38

Figure 19: LMPC control-effort analysis during closed-loop RTO–LMPC operation. The first panel compares the applied heat input 𝑄 with the RTO steady-state target 𝑄 𝑠 𝑝 . The second panel shows the LMPC correction 𝑢 0 = 𝑄 − 𝑄 𝑠 𝑝 , 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 Δ𝑄 max ∫= 5 × 103 cal/s shown as a reference. The final panel shows the cumulative absolute LMPC control effort, |𝑢 0 | 𝑑𝑡 . aggressive price variations produced larger transient deviations: the sinusoidal case resulted in a temperature RMSE of 0.761 K, whereas the random and spike profiles increased the RMSE to 1.85 K and 1.58 K, respectively. Table 11: Closed-loop performance under alternative heat-price profiles. Profile

𝐶 𝐴 RMSE

𝑇 RMSE (K)

Max. |𝑇 − 𝑇 𝑠 𝑝 | (K)

Max. |𝑢 0 | (cal/s)

Integrated positive economic gap

Max. plant violation

Baseline Sinusoidal Random Spike Ramp

6.35 × 10−4 1.10 × 10−3 3.87 × 10−3 3.41 × 10−3 7.90 × 10−4

0.278 0.761 1.845 1.580 0.299

1.972 6.345 20.216 22.053 1.880

13.81 45.66 139.31 163.46 13.16

0.702 1.451 4.238 2.149 0.000

0.00 0.00 0.00 0.00 0.00

The largest transient errors occurred for the spike and random profiles, for which the maximum temperature tracking errors reached approximately 22.05 K and 20.22 K, respectively. These cases also required larger LMPC corrections, with maximum |𝑢 0 | values of 163.46 and 139.31 cal/s. Despite these more demanding operating scenarios, the process remained within all physical bounds. The integrated positive economic gap increased from 0.702 for the baseline case to 4.24 for the random profile, indicating the expected loss in transient economic performance when the RTO target changes more abruptly. In contrast, the smooth ramp profile exhibited essentially zero positive economic gap while maintaining low tracking error. 39

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 𝐶 𝐴 tracking error, temperature tracking error, LMPC correction 𝑢 0 , 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.

E

Multi-Start Projected Adam Optimization

After training, the KKT-HardNet surrogate is used as a differentiable model for optimizing the syngas operating conditions. The optimization variables are the feed-flow ratio, reactor operating pressure, and reactor operating temperature, 𝑧 = (𝑎, 𝑃op , 𝑇op ). These variables are restricted to the same operating ranges used to generate the training data and are normalized to z ∈ [0, 1] 3 during the optimization. For a prescribed scalarization parameter 𝛼 ∈ [0, 1] , the objective is 𝐹 = 𝛼𝑈 − (1 − 𝛼) 𝑋CO2 ,

where 𝑈 denotes the total utility requirement and 𝑋CO2 = 1 − 𝐹CO2 ,out

denotes the CO2 conversion. The parameter 𝛼 controls the trade-off between utility consumption and conversion. Small values of 𝛼 place greater emphasis on maximizing conversion, whereas larger values place greater emphasis on minimizing utility usage. 40

Because the surrogate model and the resulting scalarized objective are nonlinear, the optimization problem may contain multiple stationary points. A single gradient-based run may therefore converge to different local solutions depending on the initial operating condition. To reduce this sensitivity to initialization, a multi-start strategy is used. For each value of 𝛼, the optimization is repeated from 𝑁 𝑠 randomly sampled initial operating points, and the solution with the smallest objective value is retained. This approach does not guarantee a global optimum, but it provides a more robust search of the bounded operating domain than a single initialization. Algorithm 1: Multi-start projected Adam for syngas operating-condition optimization Problem: Trained KKT-HardNet surrogate M Parameter: 𝛼 ∈ [0, 1] , 𝜂 = 0.01, 𝛽1 = 0.9, 𝛽2 = 0.999, 𝜖 = 10−8 , 𝐾 = 1000, 𝜏 = 10−8 , 𝑁 𝑠 = 100 foreach 𝛼 ∈ {0, 0.1, . . . , 1.0} do 𝐹 ∗ ← +∞ for 𝑠 = 1, . . . , 𝑁 𝑠 do Sample 𝑧 𝑘 = (𝑎, 𝑃op , 𝑇op ) uniformly within bounds and normalize to z ∈ [0, 1] 3 Initialize m0 = 0 and v0 = 0 for 𝑘 = 1, . . . , 𝐾 do Evaluate (𝐹CO2 ,out , 𝑈, . . .) = M(𝑧 𝑘 ) Compute 𝑋CO2 = 1 − 𝐹CO2 ,out and 𝐹 = 𝛼𝑈 − (1 − 𝛼) 𝑋CO2 Compute g 𝑘 = ∇z 𝐹 Update m 𝑘 = 𝛽1 m 𝑘−1 + (1 − 𝛽1 )g 𝑘 and v 𝑘 = 𝛽2 v 𝑘−1 + (1 − 𝛽2 )g2𝑘 b 𝑘 = m 𝑘 /(1 − 𝛽1𝑘 ) and b Compute m v 𝑘 = v 𝑘 /(1 − 𝛽2𝑘 ) √︁ b 𝑘 /( b Update e z 𝑘+1 = z 𝑘 − 𝜂 m v 𝑘 + 𝜖) Project z 𝑘+1 = clip(e z 𝑘+1 , 0, 1) if ∥z 𝑘+1 − z 𝑘 ∥ ∞ < 𝜏 then break if 𝐹 < 𝐹 ∗ then 𝐹∗ ← 𝐹 ∗ , 𝑇 ∗ , 𝑈∗, 𝑋 ∗ ) Store (𝑎 ∗ , 𝑃op op CO2 return Best operating conditions and utility–conversion values for each 𝛼 Each local optimization is performed using Adam. At every iteration, the KKT-HardNet surrogate is evaluated to obtain the process outputs, the scalarized objective is computed, and its gradient with respect to the operating variables is obtained. An Adam update is then performed, followed by projection onto the normalized box [0, 1] 3 . This projection guarantees that every iterate remains within the prescribed operating bounds and, therefore, within the domain represented by the surrogate training data. The optimization terminates when the change in the normalized operating variables is smaller than the prescribed tolerance or when the maximum number of iterations is reached. The complete procedure is summarized in Algorithm 1. The multi-start procedure is repeated independently for each value of 𝛼, producing a family of operating points that characterizes the trade-off between CO2 conversion and utility consumption. The projection step enforces the operating-variable bounds throughout the optimization, while the multiple random initializations reduce the dependence of the final solution on a particular starting point.

41

Related documents

• A Tutorial on Physics Constrained Machine Learning for Hybrid Modeling, Optimization, and Control#1448429
• Official Review · ICLR 2026#787029
• Comparing Classical and Quantum Machine Learning for Regression in High Energy Physics Collision Data#784846
• Comparing Classical and Quantum Machine Learning for Regression in High Energy Physics Collision Data#783256
• Safe Deep Reinforcement Learning for Building Heating Control and Demand-side Flexibility#31320