Demand Response for Flat Nonlinear MIMO Processes using Dynamic Ramping Constraints

Volatile electricity prices make demand response (DR) attractive for processes that can modulate their production rate. However, if nonlinear dynamic processes must be scheduled simultaneously with their local multi-energy system, the resulting scheduling optimization problems often cannot be solved in real time. For single-input single-output processes, the problem can be simplified without sacrificing feasibility by dynamic ramping constraints that define a derivative of the production rate as the ramping degree of freedom. In this work, we extend dynamic ramping constraints to flat multi-input multi-output processes by a coordinate transformation that gives the true nonlinear ramping limits. Approximating these ramping limits by piecewise affine functions gives a mixed-integer linear formulation that guarantees feasible operation. As a case study, dynamic ramping constraints are derived for a heated reactor-separator process that is subsequently scheduled simultaneously with its multi-energy system. The dynamic ramping formulation bridges the gap between rigorous process models and simplified process representations for real-time scheduling.


Introduction
Reducing greenhouse gas emissions requires increased renewable electricity production that, however, gives a fluctuating supply. This fluctuating supply can be compensated by consumers that react to time-varying electricity prices by shifting their demand in time in a so-called demand response (DR) (Zhang and Grossmann, 2016). DR can be attractive for energy-intensive production processes with the flexibility to modulate their production rate. To exploit the DR potential, a scheduling optimization is needed, which typically determines operational set-points for a time horizon in the order of one day (Baldea and Harjunkoski, 2014). However, such a scheduling optimization is computationally challenging for nonlinear processes that exhibit scheduling-relevant dynamics. The scheduling optimization becomes even more difficult if processes do not only consume electricity but also heating or cooling as these processes need to be scheduled simultaneously with the local multi-energy supply system (often also referred to as utility system) (Leenders et al., 2019). Local multi-energy supply systems bring integer on/off decisions into the scheduling optimization, especially as they often consist of multiple redundant units (Voll et al., 2013). Thus, the simultaneous scheduling optimization problem usually is a nonlinear mixed-integer dynamic optimization (MIDO) problem (Baader et al., 2022b).
Traditionally, such challenging scheduling MIDO problems are simplified by introducing static ramping constraints that define the rate of change of the production rate ρ as degree of freedom ν and limit ν with constant ramping limits ν min , ν max , see e.g., Carrion and Arroyo (2006); Mitra et al. (2012); Adamson et al. (2017); Zhou et al. (2017); Hoffmann et al. (2021): If additionally, the nonlinear energy demand of the production process is approximated as a piecewiseaffine (PWA) function, the complete problem can be reformulated as a mixed-integer linear program (MILP) (Schäfer et al., 2020). Traditional ramping constraints, however, have two shortcomings: They are restricted to (i) firstorder dynamics and (ii) constant ramping limits. To tackle both shortcomings, we proposed high-order dynamic ramping constraints in our previous publication (Baader et al., 2021). These high-order dynamic ramping constraints limit the δ-th derivative of the production rate ρ and use limits ν min , ν max which are functions of the production rate and its time derivatives: ρ (δ) = ν with ν min ρ,ρ, ..., ρ (δ−1) ≤ ν ≤ ν max ρ,ρ, ..., ρ (δ−1) Dynamic ramping constraints allow to represent energy-intensive processes better than static ones due to their ability to account for high-order dynamics and non-constant ramping limits that depend on the process state. In Baader et al. (2021), we demonstrated that dynamic ramping constraints can achieve solutions of the scheduling optimization problem that are close to the solutions of the original nonlinear MIDO problem while allowing to formulate MILPs that can be solved sufficiently fast for real-time scheduling. Moreover, we presented a method to derive dynamic ramping constraints rigorously for the special case of single-input single-output (SISO) processes that are exact input-state linearizable. However, this case is quite restrictive, and the question remained how to derive dynamic ramping constraints for more general processes. In particular, it was open how to apply the dynamic ramping method to multi-input multi-output (MIMO) processes that have a variable production rate ρ while, at the same time, other output variables y such as temperatures and concentrations need to be controlled using a set of control inputs u (Figure 1). The present publication extends the dynamic ramping approach to MIMO processes and presents a method to derive dynamic ramping constraints rigorously for differentially flat MIMO processes. In simple words, a MIMO process with m inputs is differentially flat if an output vector ξ of the same size m exists such that all process states and inputs can be given as a function of the output ξ and its β time derivativesξ, ..., ξ (β) (Fliess et al., 1995;Rothfuß, 1997;Rothfuß et al., 1996). We make use of the fact that a flat nonlinear model can be transformed to a linear model (Fliess et al., 1995). However, constraints that are linear for the original model, e.g., bounds on inputs, are nonlinear for the transformed model. In other words, a nonlinear model with linear constraints is transformed into a linear model with nonlinear constraints. Approximating the nonlinear bounds to the safe side with piecewise-affine functions production rate ρ control inputs u controlled outputs y t energy demand Fig. 1. A MIMO process that can vary its production rate ρ and thus its energy demand while additional control inputs u are available to control further process outputs variables y.
allows us to formulate a MILP. Note that for SISO processes, flatness is equivalent to exact input-state linearizability (Adamy, 2014), which was the main assumption in our previous work (Baader et al., 2021). For MIMO processes, exact input-state linearizability with static state transformation is a special case of flatness (Adamy, 2014).
The remaining paper is structured as follows: In Section 2, the considered demand response optimization problem and the dynamic ramping constraints are introduced. In Section 3, a method is presented to derive dynamic ramping constraints for flat MIMO processes. In Section 4, a case study is presented. Section 5 provides further discussion and Section 6 concludes the work.

Dynamic mixed-integer linear scheduling with ramping constraints
This section briefly introduces the simultaneous scheduling optimization problem of flexible processes represented by dynamic ramping constraints and multi-energy supply systems. The optimization problem (P) reads: The objective is to minimize the cumulative energy costs Φ energy at final time t f . In the following, we focus on the dynamic ramping constraints (Pb) and the process energy demand model (Pc). The remaining constraints (Pd) to (Ph) are standard constraints (Schäfer et al., 2020;. In our previous publication (Baader et al., 2021), we only considered a single input. We, therefore, only had to constrain a single time derivative of the production rate ρ with order δ, which we defined as the ramping degree of freedom ν = ρ (δ) . In the present contribution, we consider multiple inputs u and can potentially have constraints on all considered derivatives of the production rate ρ (γ) with γ = 1, ..., δ and the integer δ being the order of the highest time derivative that is constrained by input bounds. For instance, the bounds of input u 1 could directly limit the first derivative of the production rate,ρ, whereas the bounds on input u 2 could limit the second derivative of the production rate, ρ (2) , directly and only influence the first derivative,ρ, through the integration. The generalized dynamic ramping constraints (DRCs) developed in this publication therefore read: . . .
Still, only the highest considered time derivative ρ δ is a degree of freedom and thus defined as the ramping degree of freedom ν.
The limits ρ (γ) min , ρ (γ) max as well as ν min , ν max are in general nonlinear functions that we derive by coordinate transformation (Section 3). To incorporate the dynamic ramping constraint into an MILP formulation, the true nonlinear limits are approximated by piecewise-affine (PWA) functions for all considered derivatives γ = 1, ..., δ because PWA functions allow us to formulate a mixed-integer linear scheduling problem. These piecewise-affine functions need to be conservative such that they prohibit constraint violation with respect to the true nonlinear limits to guarantee that the chosen trajectory for ρ satisfies all bounds on inputs and states. Accordingly, the approximations of ρ (γ) min and ν min must always by greater than or equal to the true nonlinear limits and the approximations of ρ (γ) max and ν max must always be lower than or equal to the real nonlinear limits. Choosing the quality of the approximations allows us to explicitly balance the achievable flexibility range against the computational burden. Finally, time is discretized using collocation on finite elements to convert the dynamic ramping constraints (DRC) into a set of piecewise-affine constraints (Biegler, 2010).
The process energy demand Q dem,e (cf. Equation Pc) for an energy form e is modeled as a function of the production rate and its time derivatives: Q process dem,e ρ,ρ, ..., ρ (δ−1) , ν Similar to the DRC, a piecewise-affine function is chosen for Q process dem,e to achieve an MILP formulation. The problem formulation (P) with DRCs captures a larger share of the dynamic flexibility range of processes than static ramping constraints can do while still allowing to formulate a mixed-integer linear program. However, to choose suitable piecewise-affine ramping limits, the true nonlinear limits of the process need to be derived or approximated. In the following section, we show how these limits can be derived rigorously for MIMO processes that are differentially flat.

Deriving dynamic ramping constraints
In Section 3.1, necessary assumptions are stated and then our approach is presented in Section 3.2.

Assumptions
1. The process degrees of freedom are divided into a control input vector u and the variable production rate ρ. The process model is a system of ordinary differential equations given by: with state vector x ∈ R n , and a nonlinear right-hand side function f (x, u, ρ). We assume that there are no further inputs or disturbances to the process.
2. The control input vector u, and the production rate ρ are bounded by minimum and maximum values u min , u max , and ρ min , ρ max , respectively. Similarly, the states x have to be maintained within bounds x min , x max .
3. The process (4) is flat. That is, the process has one or multiple flat output vectors ξ. An output vector ξ is flat if it satisfies three conditions (Fliess et al., 1995;Rothfuß, 1997;Rothfuß et al., 1996): (a) The flat output vector can be given as a function φ of states x, inputs u, production rate ρ, and time derivatives of u and ρ: u,u, ..., u (α) , ρ,ρ, ..., ρ (κ) with finite integer numbers α, κ. The function φ can be seen as a transformation from the original state and input space to the flat output space. Often, it is possible to choose flat outputs that have a physical meaning, e.g., the conversion of a reactor, and that are a function of the states x only (Adamy, 2014).
(b) A backtransformation from the flat output and its derivatives to the original states x and inputs u can be found. Accordingly, the system states x and inputs u can be given as functions ψ 1 , ψ 2 of the flat outputs ξ, the production rate ρ, and a number of time derivatives of ξ and ρ: with finite integer numbers β, ζ.
Note: To check conditions (5a) -(5c), a candidate for a flat output vector ξ is needed. We assume that such a candidate for a flat output vector ξ can be identified based on engineering intuition.
4. The trajectory of the production rate ρ is determined by the scheduling optimization. Subsequently, the control input vector u ∈ R m is calculated by an underlying process control.
The flatness-based coordinate transformation is visualized in Figure 2. A nonlinear model with linear constraints is transformed into a linear model with nonlinear constraints. The number of control degrees of freedom is maintained: While, in the original model, the m inputs u k (k = 1, ..., m) are the degrees of freedom, in the linear model, the m highest time derivatives of the output components ξ (βk) k (k = 1, ..., m) are the degrees of freedom (Adamy, 2014;Fliess et al., 1995). The number of time derivatives β k can deviate between the flat output components ξ k (Adamy, 2014;Fliess et al., 1995). In the linear model, the transformed state vector Ξ is formed by the outputs and their time derivatives (except for the highest time derivative): The dimension of the transformed state vector Ξ is greater or equal to the dimension of the original state vector x (Adamy, 2014;Fliess et al., 1995). Consequently, the original nonlinear process model is converted to a linear model consisting of m integrator chains. In the linear model, every flat output ξ k can be varied with a β k -th order dynamic independently of the other flat outputs. Note that flatness is a sufficient condition for controllability of a nonlinear process (Adamy, 2014). To ease notation in the following, we introduce the ramping state vector ϕ = ρ,ρ, . . . , ρ (δ−2) , ρ (δ−1) T and its time derivativeφ

Approach
The steps to derive dynamic ramping constraints further discussed in the following subsections are summarized in Figure 3: First, a candidate for a flat output vector is selected based on two necessary flatness conditions (Section 3.2.1). Assuming this candidate is, in fact, a flat output, the original nonlinear MIMO process model is transformed into a linear MIMO model (center of Figure 3). In this linear model, the m components of the flat output vector ξ k are decoupled such that the outputs ξ k can be controlled independently of each other by manipulating the degrees of freedom ξ (βk) k ( Figure 2). However, a SISO model is needed for the dynamic ramping constraints where the production rate ρ is controlled by manipulating the ramping degree of freedom ν. Thus, second, the flat output components ξ k are coupled by choosing an operating strategy π(ρ) that defines every ξ k as function of the production rate ρ (Section 3.2.2). This coupling reduces the number of degrees of freedom from m to one, leading to a significant model order reduction. Third, backtransformations are found from the ramping SISO model to the original nonlinear MIMO model, and thus flatness is proven (Section 3.2.3).
Based on the backtransformations and the ramping SISO model, the true nonlinear ramping limits are derived from the bounds on inputs u and states x. After approximating the true nonlinear limits with piecewise-affine functions, the dynamic ramping constraints are complete and the problem (P) can be solved as MILP. 3 Find backtransformations x = ψ 1 (φ), u = ψ 2 (φ, ν) and reformulate to dynamic ramping constraints (DRCs) (Section 3.2.3) For the nonlinear MIMO model and the linear MIMO model, all variable symbols are omitted for clarity as they are identical to Figure 2. The variables of the ramping SISO model are the production rate ρ, its derivativesρ, ..., ρ (δ−1) , and the ramping degree of freedom ν.

Selection of flat output candidate and necessary flatness conditions
First, we propose to apply a necessary condition for flatness from literature based on a graph representation of the process (Schulze and Schenkendorf, 2020) to identify a candidate for the flat output vector ξ as φ x, u,u, ..., u (α) , ρ,ρ, ..., ρ (κ) (compare to condition 5a). As the second necessary condition for flatness, we propose setting up an equation system that implicitly defines the backtransformation ψ (compare to condition 5b) and check if this backtransformation equation system is structurally solvable.
For the graph representation, the process model is represented as a directed graph in which all states x i and inputs u k are represented as vertices v xi and v uk , respectively (Schulze and Schenkendorf, 2020). If ∂fi(x,u) ∂xj is non-zero, there is an edge from vertex v xj to vertex v xi (Schulze and Schenkendorf, 2020). In other words: If a state x j acts on the derivative of another state x i , an edge is drawn from x j to x i . Similarly, an edge from input u k to state x i exists, if ∂fi(x,u) ∂uk is non-zero. The necessary condition for a flat output vector ξ is that it must be possible to match the m output components ξ i to the input components u i such that there are m input-output pairs with pairwise disjoint paths through the graph that cover all process states (Schulze and Schenkendorf, 2020). An illustrative example for this necessary condition is given in the Supplementary Information (SI).
As a starting point, the states which are typically controlled can be tested as flat output candidates. Typical states to be controlled are outlet stream compositions, final effluent temperature, or the process hold-up (Jogwar et al., 2009).
Once a flat output vector candidate is identified based on the graph representation as ξ = φ(x, u,u, ..., u (α) , ρ,ρ, ..., ρ (κ) ), condition 4a is fulfilled. Further, the rank criterion needed for condition 4c can be checked easily. As a final step to show flatness, existence of the backtransformations ψ 1 , ψ 2 that give states x and inputs u as functions of the flat output ξ and its time derivativesξ, ..., ξ (β) needs to be shown (condition 4b). To obtain these backtransformations, a nonlinear system of equations needs to be developed with ξ and its derivativesξ, ..., ξ (β) on the right-hand side and left-hand side functions including the wanted quantities: states x, inputs u, and potentially derivatives of inputs. To set up the backtransformation, we make use of the fact that ξ is given by φ x, u,u, ..., u (α) , ρ,ρ, ..., ρ (κ) and the derivatives of ξ can be given as functions of x, inputs u, and derivatives of inputs by means of the total differential. The equation system for the backtransformation needs to be square such that there are as many equations as there are unknown states, inputs, and derivatives of inputs. We propose to check if this nonlinear equation system is structurally solvable by conducting an analysis similar to the structural index analysis for differential-algebraic equation systems (Unger et al., 1995) commonly used in process systems engineering. An example is given in the SI.

Operating strategy
Under the assumption that the flat output candidate identified as discussed in Section 3.2.1 is, in fact, a flat output, the nonlinear model can be transformed into a linear model. This linear model is a MIMO model with m degrees of freedom ξ (βi) i and m outputs ξ i ( Figure 2). As discussed above, the linear MIMO is transferred to a SISO ramping model by coupling the flat output components. To this end, we insert an operating strategy that gives the value of every flat output as a function of the production rate ρ. This operating strategy differentiates between two possible types of flat output components: First, flat output components might have specifications that should be maintained constant, such as outlet stream compositions, final effluent temperature, and the process hold-up (Jogwar et al., 2009). Accordingly, the operating strategy is to hold such an output ξ k constant at its nominal value such that ξ k = ξ nom k holds and thus all time derivatives are zero.
Second, flat output components may be unspecified. For instance, in our case study, one flat output component is a concentration for which no specifications are given. As every flat output component ξ k corresponds to one control input, i.e., degree of freedom, u k , if specifications are given for l outputs, and l is smaller than m, m − l flat output components remain as degrees of freedom in steady-state. Thus, the m − l free flat output component ξ k can, in principle, be chosen to have any value in steady-state as long as no variable bounds are violated. To have the optimal steady-state operating points, we use a steady-state optimization to determine the optimal values ξ k in advance as a function of the production rate ρ such that ξ k = π k (ρ). For instance, the objective can be to find the steady-state operating points that minimize energy consumption. In our case study, we choose the flat output component such that the overall heat demand is minimal for steady-state points (Section 4).
The operating strategy π k (ρ) can be any nonlinear function. The only requirement is that the function π k (ρ) must be differentiable with respect to ρ sufficiently often so that all derivatives of ξ k which are part of the backtransformation discussed in the previous section are defined by the total differential, e.g., 2 . When all output components with specifications are maintained at their nominal values and a function π k (ρ) is chosen for all other output components, the operating strategy can be summarized as ξ = π(ρ). Consequently, the flat output vector ξ and all relevant time derivatives ξ,ξ, ..., ξ (β) (compare to Equation (5b)) are defined as function of the production rate ρ and a number of its time derivatives. The highest time derivative ρ (δ) that occurs defines the order of the dynamic ramping constraint (DRC) and the ramping degree of freedom ν = ρ (δ) . Consequently, the backtransformations ψ 1 , ψ 2 discussed in the following only depend on the ramping state vector ϕ and the ramping degree of freedom ν.

Backtransformation and reformulation to dynamic ramping constraints
First, the flat output vector ξ and its time derivatives ξ,ξ, ..., ξ (β) are replaced in the nonlinear equation system for the backtransformation derived in Section 3.2.1. While ξ is replaced by π(ρ), the derivatives are replaced by building the total differential of π(ρ). Next, we solve the equation system to get x = ψ 1 (ϕ) and u = ψ 2 (ϕ, ν). It is favorable to solve the equation system analytically. Still, it is not necessary to derive the functions ψ 1 , ψ 2 analytically as they can also be evaluated numerically as long as their solution is unique. In this paper, we use the computer algebra package SymPy (Meurer et al., 2017) to obtain analytic functions ψ 1 , ψ 2 . Note that it might not be possible to solve the system of equations because the graphical and structural criteria proposed in Section 3.2.1 are only necessary flatness conditions. In case the nonlinear system of equations cannot be solved, one can test other flat output candidates. If the functions ψ 1 , ψ 2 are found, condition 3b is fulfilled, and flatness is shown at least locally on the subspace defined by the operating strategy.
This flat system constitutes one integrator chain with the degree of freedom ν that can be chosen arbitrarily at any point in time. Thus, mathematically, the production rate can be changed infinitely fast by choosing infinitely high values for ν. However, the real process control inputs u are bounded by maximum and minimum values (assumption 3), and therefore, dynamic ramping constraints (DRCs) are needed to ensure that the real process inputs are maintained within bounds.  . Constraints for ramping degree of freedom ν as function of production rate ρ for an illustrative case with first-order ramping constraints and two limiting inputs u 1 , u 2 . For first-order ramping constraints, the ramping state vector ϕ is of dimension one and equal to the production rate ρ. Consequently, the limits on the ramping degree of freedom ν only depend on ρ. The true nonlinear limits caused by the minimum and maximum values of the two inputs u 1 , u 2 are compared to static limits (dotted), linear limits (dashed) and piecewise-affine (PWA) limits (dashed-dotted).
To get the dynamic ramping constraints (DRCs), we consider the input bounds row by row. For the k-th row, we get Bounds on a state x k give an equation which is analog to Equation (9). Thus, constraints on states can be treated in the same way as constraints on inputs.
Finally, the nonlinear ramping limits need to be approximated by piecewise affine (PWA) functions. This approximation allows to explicitly balance the quality of the dynamic ramping constraints against the computational burden. In contrast to previous work (Baader et al., 2021), there can more be several inputs and states limiting the same derivative of the production rate ρ for the MIMO case. Figure 4 shows an illustrative case with first-order ramping and two limiting inputs u 1 , u 2 . The upper ramping limit is only determined by the maximum input u max 1 as the upper limit in the flat system resulting from u min 2 is always above the limit from u max 1 . In contrast, the lower ramping limit is given by an intersection between the lower limits from u min 1 and u max 2 , respectively. If static ramping limits are chosen, a large amount of the feasible region needs to be cut off for the illustrative case in Figure 4 (horizontal dotted lines). In contrast, linear (dashed lines in Figure 4) and piecewise affine (dashed-dotted lines in Figure 4) limits allow to come closer to the true nonlinear limits and thus realize a larger flexibility range. For the lower ramping limit, piecewise affine limits can be realized without the addition of binary variables as the feasible region is convex. However, for the upper limit, the feasible region is non-convex and binary variables are needed, making the optimization more computationally challenging.
4 Case study: Reactor-separator process with recycle Case study of reactor-separator process with recycle: States x are the concentrations of component A and B, c A , and c B , respectively, and the temperature T in the reactor (1) and in the flash (2). Manipulated control inputs u are the bottom stream F B , the purge stream F P , the heat input to the reactor Q 1 , and the heat input to the flash Q 2 . The scheduling degree of freedom is the production rate ρ. All other material flow rates are given as functions of F B , F p , and ρ, e.g., the reactant stream is equal to the purge F p plus the production rate ρ as no accumulation of material occurs.
In this case study, we consider a reactor-separator process with recycle consisting of a continuous stirred tank reactor (CSTR) and a flash ( Figure 5). The production rate ρ can be varied around its nominal value ρ nom between ρ min = 0.8ρ nom and ρ max = 1.2ρ nom as long as the nominal production is reached on average over the considered time horizon. A raw material A reacts to the desired product B, which can further react to an undesired product C. The process has 6 differential states: the concentration of A, c A1 , the concentration of B, c B1 , and the temperature T 1 in the reactor (1) and the analog quantities c A2 , c B2 , T 2 , in the flash (2). Apart from the production rate ρ, there are four manipulated variables: the bottom stream F B , the purge stream F p , the heat flow to the reactor Q 1 , and the heat flow to the flash Q 2 . The process equations are modified from the textbook example by Christofides et al. (2011) where 2 CSTRs and a flash are considered . Though, also the original version (Christofides et al., 2011) fulfills our assumptions, we decided to modify the example to reduce the number of states from 9 to 6 to improve readability and clarity.
The model equations comprise component and energy balances given in the SI.

Selection of flat output candidate and necessary flatness conditions
A flat output candidate ξ must have 4 components as there are four control inputs. We first consider the three states of the flash c A2 , c B2 , and T 2 , as they determine the outlet stream, and we assume that specifications for the outlet stream are given. As fourth output component, we choose the concentration c A1 based on the graph representation in Figure 6. The four input-output pairs fulfill the necessary condition for a flat output ( Figure 6).  Fig. 6. Graph representation for the reactor-separator process with recycle (compare to Figure 5). The output ξ = (c A2 , c B2 , T 2 , c A1 ) T fulfills the necessary condition for a flat output.
By differentiating the components of the output ξ = (c A2 , c B2 , T 2 , c A1 ) up to three times, we receive a structurally solvable system of equations (Table S3 in the SI).

Operating strategy
For the operating strategy, we assume that the composition and temperature of the product stream ρ must be maintained constant. Accordingly, ξ 1 = c A2 , ξ 2 = c B2 , ξ 3 = T 2 must be maintained at their nominal values ξ nom 1 = 0.4539, ξ nom 2 = 0.4610, ξ nom 3 = 455K. Thereby, the considered operating region is already significantly reduced.
As there are 4 control inputs, we can maintain the first three flat output components at their nominal values and still have one degree of freedom left. Consequently, an operating strategy for the fourth flat output ξ 4 = c A1 can be chosen freely. In a steady-state optimization, we search for the steady-state operating points that minimize the total heating Q 1 + Q 2 and fix ξ 4 to be a function of the production rate ξ 4 = π 4 (ρ). Further details on this steady-state optimization are given in the SI.

Reformulation to dynamic ramping constraints
After inserting the operating strategy defined in Section 4.2, we solve the nonlinear backtransformation equation system. To this end, we use the computer algebra package SymPy (Meurer et al., 2017) to find explicit algebraic expressions for all states and inputs, except for the temperature T 1 . The temperature T 1 has to be determined numerically because the two reaction terms in the differential Equation (S9) lead to an equation of the form Tab. 1. Structural dependency resulting from nonlinear transformation and operating strategy of states and inputs on the production rate ρ and its time derivativesρ and ν (= ρ (2) ). States which are held constant are marked with a star (*).
The implicit Equation (11) has a unique solution for the temperature, as the exponential functions are monotonic. The numerical solution of Equation (11) is found using the python package SciPy (Virtanen et al., 2020). Table 1 provides the structural dependency of the states and inputs on the production rate ρ and its time derivativesρ and ν (= ρ (2) ). The highest time derivative of the production rate that appears is two. Therefore, the ramping degree of freedom ν is equal to ρ (2) . More detailed information about the reformulations is given in the SI.
For the dynamic ramping constraints, it has to be analyzed which state and input bounds limit the time derivatives of the production rate ρ. States and inputs that are either held constant or exclusively depend on the production rate ρ do not need to be checked because it is already known from the steady-state optimization that these variables take feasible values for all considered production rates. Consequently, the variables c A,1 , c B,1 , c A,2 , c B,2 , T 2 , F B do not influence the dynamic ramping constraints (compare to Table 1).
The variables T 1 , F p , Q 2 depend on ρ andρ but not on the second time derivative ν. Accordingly, the bounds of these variables limit the first time derivativeρ. In Figure 7, the limits onρ resulting from variable bounds are shown over the production rate ρ. The calculation of these limits is explained in the SI. While the lower limit onρ results from the bound F min p , the upper limit is given by F max p for small production rates and by Q min 2 for large production rates. Graphically, we choose conservative linear functions for the limitsρ min (ρ),ρ max (ρ) (compare to (DRCb)). The heating input Q 1 is the only variable that depends on the ramping degree of freedom ν. We choose linear limits ν min (ρ,ρ), ν max (ρ,ρ) (compare to DRCd) which cover 80 % of the feasible area. Further details are given in the SI.
With the limits on ν, the second-order dynamic ramping constraints are completely parameterized and have the form:

Investigation 1: Ramp optimizations
To illustrate the ramping behavior with the derived dynamic ramping constraints, we perform two asfast-as-possible ramp optimizations shown in Figure 8. The ramp-up from minimum production rate to maximum production rate takes 1.3 hours, and the corresponding ramp-down takes 3.1 hours. The ramp-up is first constrained by the acceleration, i.e., the bounds on ν, and then by the speed, i.e., the bounds onρ. In contrary, the ramp-down is always constrained by the bounds on ν.  To visualize such an optimized ramp on the full-order process model, we simulate the ramp-up by using the optimized production rate trajectory (left part of Figure 8) as input to a simulation and calculate the control inputs u using the backtransformation function u = ψ 2 (ρ,ρ, ν) derived above. While the first three flat output components c A2 , c B2 , T 2 are maintained at their nominal values, the fourth output component c A1 follows the function of the production rate c A1 = π 4 (ρ) specified in the operating strategy ( Figure 9). All other states and inputs are within their respective bounds. Moreover, in the first half hour, when the ramp-up is limited by the limit on the acceleration ν max (compare to Figure 8, left), Q 1 is close to its maximum value (Figure 9) because the upper limit of ν max is derived from Q max 1 . However, Q 1 does not reach its maximum value due to the conservative linear approximation. During the second half hour, the ramp-up is limited by the maximum speedρ max (compare to Figure 8, left)   control input Q 2 , which limits the speed for high production rates ρ (compare to Figure 7), is close to its bounds ( Figure 9). Here, Q 2 comes very close to its bound as the conservative approximation of the ramping limit onρ is very close to the true nonlinear limit (Figure 7). Ramp optimization and corresponding simulation show that even though the dynamic ramping constraints are formed by purely linear equations, they can capture dynamics that are significantly more complex than traditional static first-order ramps. Moreover, the original process model with 6 states and 5 degrees of freedom is reduced to a dynamic ramping constraint with only 2 states, i.e., ρ,ρ, and one degree of freedom ν. Accordingly, coupling the flat outputs to the production rate reduces the model size and thus simplifies optimization.

Investigation 2: Demand response application
To demonstrate the dynamic ramping constraints in a DR application, the flexible process is considered together with a multi-energy system and additional non-flexible heat and electricity demands ( Figure  10). We consider an energy system based on . However, instead of one combined heat and power plant (CHP) and one boiler (B), we extended the system to 4 CHPs and 2 boilers to study a larger energy system with more discrete on/off-decisions. Additionally, electricity can be bought from and sold to the grid for the day-ahead price that may change hourly. More detailed information about the multi-energy system and the resulting optimization problem are given in the SI.
First, the energy system operation is optimized without accounting for the energy demand of the flexible process to obtain the energy costs of the inflexible demands only. Second, the operation of the energy system is optimized with the flexible process operating in steady-state such that the demand of the flexible process is constant (compare to Figure 11, right). Third, the DR optimization is performed using dynamic ramping constraints (compare to Figure 11, left).
The DR operation of the process reduces the total energy costs by 4.7 % compared to steady-state operation. Considering only the energy costs associated with the flexible process, the cost reduction through demand response is 10.2 %.  In the resulting operation, two types of periods can be distinguished: In times of low electricity prices, heat is preferably produced by the boilers, and electricity is bought from the grid (hours 14 -18). In times of high electricity prices, heat is preferable provided by the CHPs, and excess electricity is sold to the grid (hours 1-13, 19-24). The DR case reduces costs due to two reasons: First, the boilers are operated less. Instead of 12 hours, the boilers are only active in 9 hours ( Figure 11). The amount of heat provided by the boilers is reduced by 8 %. Second, the heat demand of the flexible process and, therefore, the electricity production of the CHPs is shifted from hours of low prices to hours with higher prices (Figures 11). For instance, the heat demand is lower in hour 15 and higher in hour 5. Consequently, the derived dynamic ramping constraints allow to reduce costs substantially compared to steady-state process operation.
The optimization problem terminates after a maximum runtime of 5 minutes that we set following . The remaining optimality gap is 4.2 %. Thus, even if the on/off-status of 6 energy system units has to be optimized simultaneously with the process operation, the optimization provides a schedule achieving substantial cost reductions within this maximum optimization runtime of 5 minutes.

Discussion
In this section, we discuss possible adaptations and limitations of our approach.
Throughout this paper, we assume that the production rate ρ is a degree of freedom (compare to assumptions in Section 3.1). However, the method can be adapted to cases where the production rate is not a degree of freedom but a component of the flat output vector ξ k . For instance, if the flow rate of the product stream would be hydraulically driven by a filling level h fill , the production rate would be given as part of the flat output by a function of h fill . In such cases, the production rate cannot be controlled directly but only through manipulating the corresponding input u k . For instance, the filling level h fill could be controlled by manipulating the feedflow of the corresponding unit.
In our operating strategy discussed in Section 3.2.2, all flat output components without given specifications are coupled to the production rate. This coupling reduces the dimensionality of the model. Instead of m flat outputs and their derivatives, only the production rate and its derivatives are variables in the optimization problem. Consequently, there are δ states and one degree of freedom. This dimensionality reduction strongly reduces the computational complexity of the optimization problem. Still, coupling all flat outputs with the production rate might be unfavorable in cases where some flat outputs can only be changed much slower than the production rate. In such cases, the outputs without specifications could be kept as independent integrator chains in the ramping model (right part in Figure 3). In our case study, the concentration c A1 , which is the fourth flat output component, could be uncoupled such that there are two integrator chains in the ramping model: one for the production rate ρ and one for the concentration c A1 . Consequently, dynamic ramping constraints need to be derived for both integrator chains. While this uncoupled version makes the ramping model computationally more challenging, it is also more flexible and thus might enable higher profits in some cases. An extreme example is electrolyzers that can often adapt their production rate rapidly but have slow temperature dynamics (Simkoff and Baldea, 2020;Flamm et al., 2021). Thus, for electrolyzers, it is not favorable to couple the temperature with the production rate. Instead, we have shown recently in a conference paper that it is favorable to keep both production rate and temperature as degrees of freedom and formulate dynamic ramping constraints on the temperature (Baader et al., 2022a).
The present work focuses on demand response applications where the production rate ρ is the main scheduling-relevant variable. However, it is straightforward to adapt the approach to cases where a different variable is scheduling-relevant. For example, in multi-product processes, the concentration may be varied to yield different products. Thus, the scheduling needs to account for the dynamics of the concentration during product transitions (Flores-Tlacuahuac and Grossmann, 2006;Baader et al., 2022b). Our approach can be transferred to such a multi-product process if the production rate ρ is replaced by the concentration.
A limitation of this work is that solving the equation system to derive the backtransformation requires a lot of manual effort, even though a computer algebra system is used as support. It is an open question to what extent our approach can be automated to enable the analysis of large-scale processes.
Additionally, the assumption of flat process models limits the applicability of our approach. For nonflat process models, the ramping limits could still be derived as functions of all process states x. However, it is not possible to find the coordinate transformation from the process states x to a transformed state vector Ξ as defined in Equation (6), and thus, the ramping limits cannot be given as functions of the ramping state vector ϕ only. Hence, an extension to non-flat processes would be partly heuristic as the state vector x in the ramping limits would have to be approximated based on the ramping state ϕ.

Conclusion
Dynamic ramping constraints simplify the simultaneous demand response scheduling optimization of production processes and their energy systems compared to an optimization considering the full-order process model. Still, dynamic ramping constraints can capture more flexibility than traditional static ramping constraints as they allow high-order dynamics and non-constant ramp limits. In this paper, we extend our method to rigorously derive dynamic ramping constraints from input-state linearizable single-input single-output (SISO) processes to flat multi-input multi-output (MIMO) processes. In the MIMO case, dynamic ramping constraints reduce the problem dimensionality by coupling all flat outputs with the production rate. In our case study, a system with 6 states and 5 degrees of freedom is reduced to 2 states and one degree of freedom. Additionally, we demonstrate that an operational strategy can be chosen for flat outputs such that the steady-state production points are optimal, e.g., with respect to energy consumption.
Our case study demonstrates that dynamic ramping constraints allow for finding DR schedules for flexible processes and multi-energy systems that substantially reduce energy costs compared to a steady-state operation. Even though discrete on/off-decisions in the multi-energy system add to the computational complexity, the problem can be solved within the time limit for online scheduling.
Overall, dynamic ramping constraints allow bridging the gap between nonlinear process models and simplified process representations suitable for real-time scheduling optimization.

Declaration of Competing Interest
We have no conflict of interest.  In this section, the two necessary conditions for flatness (compare to Section 3.2.1 in the main paper) are illustrated using an example (E) with three states and two inputs:

Nomenclature
In the example (E), the output vector ξ = (x 1 , x 2 ) T does not fulfil the necessary graphical condition for flat outputs because the state x 3 is not covered ( Figure S1, center). However, the vector ξ = (x 3 , x 2 ) T is a candidate for a flat output as all states are covered ( Figure S1, right). To receive the equation system for the backtransformations, the flat output candidate ξ = φ (x) = x 3 x 2 has to be differentiated two times to obtain a square backtransformation equation system: Structurally, the backtransformation equation system (S1) -(S3) is solvable for the three states x 1 , x 2 , x 3 , the two inputs u 1 , u 2 , and the input derivativeu 2 (Table S1). The input derivativeu 2 is part of the equation system as both second derivatives ξ 1 , ξ (2) 2 are functions ofu 2 .

Illustrative example for multiple feasible regions
As discussed in Section 3.4 of the main manuscript, multiple distinct feasible regions can occur for the dynamic ramping constraints. To illustrate this phenomenon, we consider an illustrative model (IM) with Tab. S1. Analysis of structural solvability of the backtransformation equation system to calculate states x, and inputs u with flat output candidate ξ = (x 3 , x 2 ) T for example process (E). All states and inputs that appear in an equation are marked with a cross. A necessary condition for solvability is that it must be possible to circle exactly one cross in every row and exactly one cross in every line. In other words: Each equation should be solved for one variable.
states inputs and derivatives three states x, two inputs u, and a production rate ρ: As shown in the following, this process has a flat output ξ = φ (x) = x 1 x 3 . We assume an operating strategy ξ = π(ρ) = aρ bρ with two constants a, b.
For x 1 , and x 3 , the back transformation is given by the function π. For x 2 , Equation (IMa) is set equal toξ 1 =ẋ 1 = aρ which gives Similarly, insertingẋ 2 = aρ (2) −ρ into Equation (IMb) gives Finally, insertingξ 2 =ẋ 3 = bρ into Equation (IMc) gives As all states and inputs can be given as functions of production rate ρ and derivativesρ, ρ (2) , the process is de facto flat for the given operating strategy. Two feasible regions for the dynamic ramping constraints can result from the second input u 2 being a quadratic function of the first derivative of the production rateρ (Equation (S6)) ( Figure S2). One could either restrict the operating range to one of the two regions or would have to introduce a binary variable that indicates which region is active. 3 Supplementary information to case study

Model equations
In this section, we present the model equations of the case study (compare to Figure 4 in the main paper). All symbols other than states x = (c A1 , c B1 , T 1 , c A2 , c B2 , T 2 ) T , control inputs u = (F B , F p , Q 1 , Q 2 ) T , and production rate ρ are parameters marked in gray with values listed in Table S2. The index 0 indicates the feed quantities.ċ

Operating strategy
As discussed in the main paper, the fourth flat output ξ 4 can be chosen freely. Here, we calculate the optimal steady-state operating points for ξ 4 and fix ξ 4 to be a function of the production rate ξ 4 = π 4 (ρ). Accordingly, we minimize the total heat input to the process Q 1 + Q 2 as the objective function to obtain the steady-state operating points. First, we determine the optimal steady-state operating points.
To this end, a steady-state optimization problem (P s ) is formulated by sampling the production range ρ min ≤ ρ ≤ ρ max with 21 equally distributed points ρ j . The optimization problem (P s ) reads: For every point j, the process needs to be in steady state (Equation (P s b)), the first three flat output components ξ k = φ k (x) need to be at their nominal values (Equation (P s c)), and inputs and states needs to be within bounds (Constraints (P s d), (P s e)). All bounds are given in Table S4. We implement the optimization problem P s using pyomo (Hart et al., 2017) and solve it using BARON version 20.10.16 (Khajavirad and Sahinidis, 2018) in heuristic mode. Next, a function π 4 (ρ) needs to be chosen that does not increase the objective function (P s a) too much. In other words, the steady-state operation points reached with the additional constraint should give an objective function close to the objective function reached without the additional constraint (P s f ). As a first candidate, we investigate the simplest possible operating strategy, that is to hold ξ 4 constant, i.e., π 4 (ρ) = ξ nom 4 . This optimization leads to a feasible solution. Accordingly, it is feasible to always hold the fourth flat output component at a constant value. However, the objective function worsens by 2.9 % compared to the first optimization without constraint (P s f ) because it is not possible to choose a single constant value ξ nom 4 that is optimal for all operating points. As a second candidate, Tab. S4. State and input bounds variable lower bound upper bound we therefore investigate the linear strategy with optimization degrees of freedom a ξ4 0 , a ξ4 1 . The optimization finds values a ξ4 0 , a ξ4 1 such that the objective value worsens only 0.05 % compared to the first optimization without the coupling constraint (P s f ). As this linear strategy is very close to the theoretical optimum, we do not consider other functions, e.g., higher-order polynomials, and continue with the linear operating strategy (S13).

Derivation of flat process reformulation
In this section, the nonlinear coordinate transformation is derived based on the nonlinear system of equations and the operating strategy. The aim is to express the states x = (c A1 , c B1 , T 1 , c A2 , c B2 , T 2 ) T and control inputs u = (F B , F p , Q 1 , Q 2 ) T as functions of production rate ρ, and its first two derivativeṡ ρ, and ν.
For the reformulation, first, the strategy for the 4 flat output components ξ = (c A2 , c B2 , T 2 , c A1 ) T discussed in Section 4.2 of the main manuscript is inserted: (S15) While the first three outputs are chosen independent of the production rate, the fourth output c A1 is a function of the production rate ρ. Next, we consider the equation for the derivative of the first flat output, i.e.,ξ 1 =ċ A2 . As c A2 is held constant, the derivativeċ A2 is zero and we solve Equation (S10) for the bottom flow F B that is a function of the production rate ρ, i.e., , and c A1 (ρ) as defined in Equation (S17). The second output is maintained constant, too. Accordingly, we receive the conditionξ 2 =ċ B2 = 0. After inserting F B (ρ) from Equation (S18) into Equation (S11), we get the equation: As ρ is always nonzero, we get c B1 as: From the third output component, we get the conditionξ 3 =Ṫ 2 = 0 and thus solve Equation (S12) for Q 2 : The input Q 2 still depends on the temperature T 1 for which no expression has been derived up to this point. From the fourth flat output component, we getξ 4 =ċ A1 = a ξ4 1ρ (compare to Equation (S17)) and thus Equation (S7) changes to: Equation (S22) is solved numerically for T 1 . However, first, an expression for the input F p is needed. This expression is derived from the condition that according to our operating strategy the second derivative of the second output is zero, i.e., ξ (2) 2 = c (2) B2 = 0. From this condition, we receive with a nonlinear function Ψ Fp (ρ, T 1 ), which is not stated here to preserve readability. After inserting Equation (S23) into Equation (S22), Equation (S22) implicitly defines T 1 as a function of ρ andρ, T 1 (ρ,ρ).
As discussed in the main paper, we cannot solve for T 1 analytically, but the solution received numerically is unique. With T 1 (ρ,ρ), we can calculate Q 2 as Q 2 (ρ,ρ) (compare to Equation (S21)) and F p as F p (ρ,ρ) (compare to Equation (S23)). The only variable that is missing is the reactor heating Q 1 . To get an expression for Q 1 , Equation (S22) is differentiated with respect to time. Before doing so, we replace F p with Ψ Fp (ρ, T 1 ) (compare to Equation (S23)) and insert the functions for F B and c A,1 , which depend on ρ. After bringing everything to the left hand side, Equation (S22) changes to 0 = Ψ(ρ,ρ, T 1 ).
Building the total differential with respect to time gives and Ψ 0 (ρ,ρ, T 1 ) as the remaining part.
As Q 1 enters linearly, we can easily solve Equation (S25) for Q 1 . Moreover, the ramping degree of freedom ν enters linearly as well, which allows to solve Equation (S25) for ν and calculate the bounds of ν from the bounds of Q 1 (compare to Figure S3). At this point, all states and inputs can be calculated as functions of the production rate ρ and its first two derivativesρ, and ν.
Note that we do not derive expressions for the input derivativesḞ B , F B , andḞ p that are shown in Table S3. Instead, we replace F B with F B (ρ) (Equation (S18)), and F p with F p (ρ, T 1 ) (Equation (S23)) in the equations where they appear and subsequently differentiate the complete equations with respect to time.

Calculating the limits forρ
As discussed in the previous section, the limits for ν can be calculated in a straightforward way from Equation (S25). In this section, we calculate the limits forρ resulting from the bounds of T 1 , Q 2 , F p ( Figure 6 in the main paper). First, the bounds are calculated as a function of T 1 by dividing Equation (S22) by a ξ4 1 , givingρ asρ = θ T1 (ρ, T 1 ), with a nonlinear function θ T1 . Note that while Equation (S22) cannot be solved for T 1 analytically, Equation (S22) can be solved analytically forρ. Accordingly, inserting the bounds of T 1 gives the limits ofρ shown in Figure 6 in the main paper.
In the case of Q 2 , which is calculated as Q 2 (ρ, T 1 ) from Equation (S21), Equation (S21) cannot be solved forρ as no analytic expression was derived for T 1 that depends on ρ andρ. Instead, Equation (S21) is solved for T 1 as T 1 (Q 2 , ρ), which is straightforward as T 1 enters Equation (S21) linearly. Inserting which allows to calculate the limits ofρ from the bounds of Q 2 as shown in Figure 6 in the main paper. For the purge stream F p , we follow the same strategy to first calculate T 1 and thenρ. However, Equation (S23) cannot be solved for T 1 as function of F p and ρ analytically due to the exponential functions (see discussion around Equation (11) in main paper). Consequently, Equation (S23) is solved numerically for T 1 as function of the bounds on F p . This value of T 1 is then used to calculate the limits ofρ shown in Figure 6 in the main paper.

Ramping limits on ν
In this section, the derivation of the ramping limits on the ramping degree of freedom ν is given in detail. First, the space given by the limits ρ min , ρ max andρ min (ρ),ρ max (ρ) is sampled using 51×51 points. For every point, the true limits ν min true (ρ,ρ), ν max true (ρ,ρ) are calculated from Q min 1 , Q max 1 (compare to Equation (S25)) ( Figure S3). For a conservative approximation of these limits, we test linear functions ρ / ρ n o m with parameters a min 0 , a min ρ , a miṅ ρ , a max 0 , a max ρ , a maẋ ρ . The parameters are calculated in an optimization, which minimizes the difference to the true bounds while ensuring conservativeness such that ν min true (ρ,ρ) ≤ ν min (ρ,ρ) < ν max (ρ,ρ) ≤ ν max true (ρ,ρ) holds for every point (ρ,ρ). To evaluate the quality of the bounds, we calculate the coverage C(ρ,ρ): which is the distance between the chosen bounds in relation to the distance between the true bounds. The arithmetic mean of the coverage is 80 %. Accordingly, the chosen linear limits capture approximately 80 % of the feasible region while being computationally cheap.
3.6 Multi-energy system and optimization problem for demand response application In this section, we describe the multi-energy system and the optimization problem in the second study, which considers a demand response application. Our aim for this case study is to construct an illustrative example that demonstrates the application and capabilities of our new method. The considered multienergy system is based on the third day of the benchmark problem given by . Note that the data are publicly available at Helmholtz Energy Computing Initiative (2021). The given heat and electricity demands are multiplied by 3 as we increase the number of units from 2 to 6. Instead of one combined heat and power plant (CHP) and one boiler (B) as in , we use 4 CHPs and 2 boilers. All 4 CHPs have a nominal thermal output power of 450 kW, and the 2 boilers have 530 kW. The nominal efficiencies are given in Table S5.
Tab. S5. Nominal efficiencies of energy-system components. The minimum part-load fraction is 50% for the CHPs and 20% for the boilers . Following Voll (2014), the part-load efficiency curves are discretized with one affine element to achieve an accurate discretization.
The electricity price used is the German day-ahead price from April 2nd, 2021 (SMARD Marktdaten, 2021). This price curve is chosen because it features times with low prices that make the operation of the boiler favorable and times of high prices which make the operation of the CHPs favorable.
For the optimization problem, which schedules the flexible operation of the reactor-separator process, a process energy demand model is needed (Pc). As the dynamic ramping constraint is second order, the heat demand of the process is modeled as a function of ρ,ρ, and ν, i.e., Q process dem,heat (ρ,ρ, ν).
To fit the heat demand model, we sample the operating region given by the limits ρ min , ρ max ,ρ min (ρ), ρ max (ρ), and ν min (ρ,ρ), ν max (ρ,ρ) using 11×11×11 points. First, a purely linear function is used for the heat demand and the parameters are fitted using linear regression. Because the average absolute approximation error is 10% of the nominal heat demand, we try a piecewise-affine (PWA) model to increase the accuracy. For this new PWA model, the operating region is split into segments in which affine functions are used. Here, the operating region is split into four segments by making two distinctions: Ifρ is greater or smaller than zero and if ν is greater or smaller than zero. To implement the PWA model in the optimization, two binary variables are needed: zρ indicates whetherρ is positive and z ν indicates whether ν is positive. No distinction is made for ρ as we find that the heat demand is roughly linear with ρ. The PWA model reduces the mean approximation error to 5% of the nominal heat demand. The optimization problem is discretized using orthogonal collocation with 2 elements per hour and 3 collocation points per element. The binary variables are discretized with one-hour timestep to match the hourly changing prices and demands. There are 8 binary variables per timestep: 4 on/off binaries for the CHPs, 2 on/off binaries for the boilers, and 2 binaries for the PWA heat demand model. Consequently, the optimization problem has 8×24=192 binary variables for the 24 h time horizon. The ramping degree of freedom ν is discretized to be piecewise linear for every hour.
The optimization problem is formulated using pyomo (Hart et al., 2017(Hart et al., , 2011 and pyomo.dae (Nicholson et al., 2018) for discretization. The optimization problem is solved using gurobi version 8.1.0 (Gurobi Optimization, LLC, 2021). As  state that the maximum acceptable optimization runtime for scheduling problems is typically between 5 and 20 minutes, we set the maximum optimization runtime to 5 minutes. All calculations are performed on a Windows 10 machine with an Intel(R) Core(TM) i5-8250U CPU and 24 GB RAM.