**2. Background**

#### *2.1. Multi-Parametric Programming*

Overall, multi-parametric programming problems are concerned with the effect of parametric perturbations on the optimal solution of an optimisation problem. Consider the following optimisation problem:

$$z(\theta) = \min\_{\substack{\mathbf{x} \in \mathbb{R}^{n\_{\mathbf{x}}} \\ \text{subject to }}} f(\mathbf{x}, \theta)$$

where *θ* stands for the vector of uncertain parameters, which is *<sup>n</sup>θ*-dimensional, **x** is the *nx*-dimensional vector of continuous decision variables, **g**(**<sup>x</sup>**, *θ*) is the vector of inequality constraints and *f* is the objective function as a mapping R*nx*×*nθ* → R, both of which are assumed to be C<sup>2</sup> (twice continuously differentiable). Problem (1) is a multi-parametric program and its solution results in the partition of the parametric R*<sup>n</sup>θ* -space into a number of regions, also know as critical regions (CRs). Within each CR, the optimal solution and the objective value are given as functions of the uncertain parameters, i.e., **x**(*θ*) and **<sup>z</sup>**(*θ*), respectively. Even though mp-P has been studied quite actively, the class of multi-parametric nonlinear programming problems remains one of the most challenging ones [11,12]. Depending on the convexity of the nonlinear functions that form Problem (1), different solution techniques have been proposed in the literature to date.

Advances in the algorithms and theory of parametric nonlinear programs (p-NLPs) date back to the early works of Fiacco [13] and Bank et al. [14]. More specifically, in the books of Bank et al. [14] and Fiacco [13], a collection of the early research works for parametric NLPs can be found and invaluable theoretical foundations for some classes of convex p-NLPs with perturbations in the OFC and the right-hand side of the constraints are provided. Even though the term "parametric nonlinear optimisation/programming" was widely established from the aforementioned works, the early works on numerical stability analysis of NLPs by [15,16] and the work of Robinson [17,18] on generalised equations provided a significant way of studying the effect of parametric variations on the optimal solution of p-NLPs. Kyparisis [19] studied the uniqueness and differentiability of solutions of parametric nonlinear complementarity problems while in Ralph and Dempe [20], the

directional derivatives of parametric nonlinear programs were used to characterise their explicit solution. However, the first algorithm for the multi-parametric case of convex NLPs was due to Dua and Pistikopoulos [21]. The authors, based on the findings about the convexity properties of the parametric value function (*z*(*θ*)), devised an iterative procedure in which the integer variables were fixed by the solution of a primal mixed integer nonlinear program (MINLP) and the resulting mp-NLP was then transformed into an mp-LP following the outer approximation idea. Because of the value function's convexity property, the maximum error of the approximation occurs at the vertices of the CRs and if the error is greater than the prespecified tolerance, the CR is partitioned again; otherwise, integer and parametric cuts are implemented and then the algorithm iterates until the primal MINLP is infeasible. The same algorithm was revisited by Acevedo and Salgueiro [22], where the authors proposed heuristics to improve its computational efficiency while quadratic approximations were studied by Johansen [23] and Domínguez and Pistikopoulos [24]. An approximate algorithm for the solution of convex mp-NLPs was proposed by Johansen [25], who proposed the consecutive subdivision of the parametric space in hyper-rectangles and the interpolation of the parametric solution through the solution of 2*<sup>n</sup>θ* NLPs at each step. Further approaches involve the geometric vertex search by Narciso [26] and sub-gradient methods by Leverenz et al. [27]. For the nonconvex cases, Dua et al. [28] developed suitable parametric under/overestimators which were then incorporated into a spatial branch and bound routine for the global optimisation of the nonconvex problem within −tolerance. For a more thorough discussion on the algorithms that have been proposed for the solution of mp-NLPs, the interested reader is directed to the review of Domínguez et al. [29] while Hale [30], in her doctoral thesis, also offers a thorough discussion on several classes of parametric optimisation. Fotiou et al. [31,32] initially studied the polynomial multi-parametric programming problem with application to control, however, their approach did not include the definition of final nonconvex CRs, while the mixed integer polynomial case was studied by Charitopoulos and Dua [33] and a procedure for the computation of exact nonconvex CRs was presented. Despite the aforementioned research effort, mp-NLP problems remain one of the most difficult to tackle and, as illustrated in Table 1, all the aforementioned algorithms can handle uncertain parameters only on the right-hand side (RHS) of the constraints. Recently, Pappas et al. [34], by generalising the basic sensitivity theorem of Fiacco [13], devised an algorithm for the exact solution of multi-parametric quadratically constrained quadratic programs.


**Table 1.** Summary of multi-parametric nonlinear programmming algorithms.

Among the wide range of appplications that multi-parametric programming has been applied to, the invention of explicit model predictive control (mp-MPC) is undoubtedly the most dominant area where mp-P has had the biggest impact [10,12,35]. The main concept of mp-MPC is that instead of solving the optimisation problems related to standard MPC at each sampling instance, the state of the system is treated as an uncertain parameter and an mp-P can be solved offline to derive the explicit control solution once and for all [10,36].

The general formulation of mp-MPC for discrete time systems is shown by (2):

$$(mp-MP)\begin{cases} & \Phi(\mathbf{x}(t\_k)) = \min\_{\mathbf{u}} \sum\_{i=0}^{P\_{\text{fl}}-1} \mathcal{L}(\mathbf{x}\_{t\_i}\mathbf{u}\_t) + E(\mathbf{x}\_N) \\ & \text{subject to:} \quad \mathbf{x}\_{t|t=0} = \mathbf{x}(t\_k) \\ & \mathbf{x}\_{t+1} = f(\mathbf{x}\_t\mathbf{u}\_t) \\ & y\_{t+1} = h(\mathbf{x}\_t, \mathbf{u}\_t) \\ & \mathbf{A}\mathbf{x}\_t \le \mathbf{a} \qquad t = 0, 1, \dots, P\_H - 1 \\ & \mathbf{B}y\_t \le \mathbf{\mathcal{B}} \qquad t = 0, 1, \dots, P\_H \\ & \mathbf{C}\mathbf{u}\_t \le \mathbf{\mathcal{D}} \qquad t = 0, 1, \dots, P\_H \\ & \mathbf{C}\mathbf{u}\_t \le \mathbf{\mathcal{D}} \qquad t = 0, 1, \dots, P\_H \end{cases} \tag{2}$$

where **x***t*, **u***t*, **z***t* are the state, control input and system output vectors, respectively, at every sampling point, t, and are *nx*, *nu*, *ny*-dimensional. **A**, **B**, **C** are matrices of appropriate dimensions and *α*, *β*, *γ* vectors of pertinent dimensions which represent inequality constraints for the state, output and control inputs while L: R*nx*+*nu* → R is a stage cost and *E* : R*nx* → R is a terminal cost function. The repetitive solution of Problem (2) provides the optimal cost <sup>Φ</sup>(**x**(*tk*)) and the optimisation vector, i.e., the sequence of optimal control inputs **u**<sup>∗</sup> <sup>=</sup>**u**<sup>∗</sup>1, **<sup>u</sup>**<sup>∗</sup>2,..., **<sup>u</sup>**<sup>∗</sup>*PH*−<sup>1</sup> over the finite prediction horizon *PH*. Compared to the conventional model predictive control fashion, in which an optimisation problem is solved at each sampling point, through the mp-MPC notion, the explicit control law is calculated offline once and for all. The solution of the resulting mp-P problem results in the optimal control inputs as explicit functions of the (uncertain) parameters, i.e., the state of the system at each sampling instance, along with the corresponding CRs, as shown by Equation (3).

$$\mathbf{u}^\* = \begin{cases} \nu\_1(\mathbf{x(t\_k)}) & \text{if } \mathbf{x(t\_k)} \in \mathsf{CR}\_1 \\ \nu\_2(\mathbf{x(t\_k)}) & \text{if } \mathbf{x(t\_k)} \in \mathsf{CR}\_2 \\ \vdots & \vdots \\ \nu\_\omega(\mathbf{x(t\_k)}) & \text{if } \mathbf{x(t\_k)} \in \mathsf{CR}\_\omega \end{cases} \tag{3}$$

For the case of MPC for linear systems, instead of solving a quadratic program at each sampling instance, the explicit MPC requires the offline solution of an mp-QP while online, so only simple function evaluations are required [10,37,38]. This concept is also known as online via offline optimisation and it is shown in Figure 3.

**Figure 3.** Online via offline optimisation framework [10,39].

Even though mp-MPC is the niche area of mp-P, designing mp-MPCs of nonlinear systems for set-point tracking is still a computationally strenuous task as one has to design an mp-MPC for each set-point based on the algorithms that exist in the literature to date [40]. Next, we review two mathematical techniques that will enable the development of novel multi set-point explicit controllers though the algorithm we propose in the present work.

#### *2.2. Computer Algebra*

#### 2.2.1. Grobner Bases Theory ¨

The idea of the present work is to devise an algorithm for the solution of mp-NLPs from an analytical and not numerical perspective, and the reason is two-fold. Firstly, for the case that we are interested in, i.e., mp-NLPs with combined uncertainty on the RHS and OFC, no theoretical foundations exist for the computation of the explicit solution like the basic sensitivity theorem of Fiacco [13] which serves as the basis of the numerical mp-P approaches. Secondly, because of the nonconvex nature of the parametric problem, the numerical solution would require global optimisation techniques similar to the ones presented in Dua et al. [28] and would lead to an explosion in the number of convex approximate CRs.

Gro¨bner bases theory was introduced by the doctoral research of Bruno Buchberger [41] as a method of analytically solving systems of multivariate polynomial equations. In brief, Gro¨bner bases and the Buchberger algorithm can be considered as the polynomial counterpart of Gaussian elimination for the case of nonlinear systems. Before formally defining what a Grobner basis is, it would be useful to define some preliminary concepts. ¨

#### **Definition 1.** *Power products*

*Let* R *be any field and let* <sup>R</sup>[*<sup>x</sup>*1,..., *xn] be the ring of polynomials in n-indeterminates. Any polynomial can be described as a sum of terms of the form: αxβ*<sup>1</sup> 1 ··· *xβ<sup>n</sup> n with α* ∈ R *and βi* ∈ N, *i* = 1, . . . , *n and the term xβ*<sup>1</sup> 1··· *xβ<sup>n</sup> n is called a power product.*

#### **Definition 2.** *Term order*

*A term order is defined with regard to a set of power products (Tn) and imposes a total order < on the set in compliance with the conditions below:*


A number of alternative power product orderings exist but the most commonly employed is the the lexicographic one due to its computational efficiency [42]. Lastly, the notion of *ideals* is crucial within the Grobner bases theory. ¨

#### **Definition 3.** *Ideals*

*Let* R *be a field and* <sup>R</sup>[*<sup>x</sup>*1, *x*2,..., *xn*] *be a ring over the field of n-variate polynomials. Let a finite subset of the field, G* = {*g*1, *g*2,..., *gt*}*, then an ideal I can be generated by G as follows:*

$$I\_{\perp} = \{ \sum\_{i=1}^{n} \mu\_i \mathbb{g}\_i \mid \mu\_i \in \mathbb{R}[\mathbf{x}\_1, \mathbf{x}\_2, \dots, \mathbf{x}\_n], \mathbb{g}\_i \in \mathbb{R}, \forall i \}$$

For problems that accept analytical solutions, their existence is guaranteed by the Hilbert basis theorem [43], which also guarantees that algorithms used to compute Gro¨bner bases can terminate in a finite number of steps.

#### **Definition 4.** *Grobner basis [* ¨ *41]*

*A set of non-zero polynomials G* = {*g*1,..., *gt*}*, contained in an ideal I, is called a Grobner* ¨ *basis for I if and only if for all g* ∈ *I, such that g* = 0*, there exists i* ∈ {1, . . . , *t*} *such that l p*(*gi*) *divides l p*(*g*)*, where l p*(·) *stands for the leading power product of a polynomial function.*

For the calculation of Gro¨bner bases, apart from Buchberger's algorithm, Fauge`re has presented algorithms F4 [44] and F5 [45] as two variants of another algorithm. They exploit concepts from linear algebra and represent polynomials using matrix forms, thus enabling successive truncated Gro¨bner bases to be created. Lastly, software implementations of different algorithms that compute Gro¨bner bases computations can be found in freely available CAS such as Singular, SymPy and SageMath as well as commercial tools like Maple and Mathematica.

#### 2.2.2. Cylindrical Algebraic Decomposition (CAD)

The notion of cylindrical algebraic decomposition was presented by Collins in 1975 [46] in an effort to solve the problem of quantifier elimination over real closed fields. In this article, as will be shown later on, CADs are used for computing nonconvex regions in the space of parameters. Thus, we provide the following definitions for ease of exposition in the algorithmic steps that we detail later on in the manuscript.

#### **Definition 5.** *Semi-Algebraic Sets [43].*

Let <sup>R</sup>[*<sup>x</sup>*1, *x*2,..., *xn*] indicate the ring of polynomials in *n*-indeterminates with real coefficients. If, for example, a subset S of R*n* can be developed by a finite number of applications of the complementation, union and intersection operations, it is called semialgebraic and can have the following form:

{**x** ∈ R*n*| **g**(**x**) ≤ <sup>0</sup>}, where **g** ∈ R[*X*]

#### **Definition 6.** *Standard atomic formula*

*A formula that includes a functional relation over a polynomial ring in either of the ways shown below is called a standard atomic formula:*

$$\text{g(x)} = 0, \text{ g(x)} \neq 0, \text{ g(x)} < 0, \text{ g(x)} > 0, \text{ g(x)} \leq 0, \text{ g(x)} \geq 0$$

**Proposition 1** ([43])**.** *Semi-algebraic sets of* R*n can be written as a finite union of semi-algebraic sets of the form:*

$$\{\mathbf{x} \in \mathbb{R}^{\mathbb{N}\times} \mid \mathbf{g}\_1(\mathbf{x}) = \dots = \mathbf{g}\_{\omega}(\mathbf{x}) = 0, \mathbf{g}\_{\omega+1}(\mathbf{x}) > 0, \dots, \mathbf{g}\_l(\mathbf{x}) > 0\}$$

*where g*1 ..., *g<sup>ω</sup>*, *gω*+1,..., *gt are in g* ∈ <sup>R</sup>[*X*]*.*

The proof of the proposition can be found in the book of Bochnak et al. [43].

Using the definitions and propositions given above, in summary, one can use CAD routines to compute the solution to polynomial inequalities. In the process of computations, one partitions the related space over finite cells and qualifies whether or not standard atomic formulas hold. A comprehensive exposition on the solution of polynomial inequalities using CAD is given at the book of Jirstrand [47].

#### **3. Algorithms and MPC**

#### *3.1. Multi Set-Point Explicit Controller via Multi-Parametric Programming*

In the context of multi-parametric model predictive control, the state of the system at each sampling point is treated as an uncertain parameter and as a result an mp-P with RHS uncertainty arises [10,37,38,40]. Its solution results in the explicit control law, i.e., the control decisions as explicit functions of the system's initial conditions at a sampling instance along with the related CRs.

In many applications, however, particularly those related to continuous manufacturing, there is a grea<sup>t</sup> need for fast calculations in order to communicate decisions between the different layers of decision making in an effective manner. For instance, set-point tracking goals for APC are, most of the time, passed down from either the functionality of process scheduling or RTO [39,48]. In these cases, it becomes obvious that explicit MPC

can provide a significant advantage in computational time by treating the set-points or estimated model inputs as uncertain parameters. One way to design such explicit controllers, assuming the existence of the *nu* set-point, is to solve *nu* mp-P problems and thus design *nu* mp-MPCs. Another way, which has not been investigated in the literature, is to consider the set-points and/or model inputs as uncertain parameters and thus derive a multi set-point explicit MPC. Conceptually, by doing so, we would design a layered controller as given by Figure 4.

Conventionally, when mp-MPC is employed for set-point tracking of nonlinear systems, one would have to compute a different mp-MPC for each of those set-points as well as account for any time delay in the offline solution of the related mp-P should a new set-point arise. Stringent market regulations and an increasingly volatile market environment lead process industries to constantly optimise their operations and give rise to new set-points from a control perspective which in turn hinder the deployment of mp-MPC. As illustrated by Figure 4, in this article, by considering the set-points as uncertain parameters which lie within prespecified bounds, we overcome the abovementioned drawback of explicit MPC since, by solving one mp-NLP, we can design a "multi set-point" mp-MPC for nonlinear systems.

**Figure 4.** Concept of a multi set-point mp-MPC setting where instead of separate explicit controllers one designs a universal explicit MPC for all the set-points.

The design problem of a multi set-point mp-MPC can be formulated as in (4).

$$\begin{aligned} \mathcal{Y}(\mathbf{x}(t\_k), \mathbf{x}\_{sp}) &= \min\_{\mathbf{u}} \sum\_{t=0}^{P\_{\mathbf{f}}-1} \mathcal{L}(\mathbf{x}\_t, \boldsymbol{\mu}\_t, \mathbf{x}\_{sp}) + E(\mathbf{x}\_{N\_t}, \mathbf{x}\_{sp}) \\ \text{subject to:} \quad & \mathbf{x}\_{t|t=0} = \mathbf{x}(t\_k) \\ & \mathbf{x}\_{t+1} = f(\mathbf{x}\_t, \boldsymbol{\mu}\_t) \qquad t = 0, 1, \dots, P\_H - 1 \\ & \mathbf{y}\_{t+1} = h(\mathbf{x}\_t, \boldsymbol{\mu}\_t) \qquad t = 0, 1, \dots, P\_H - 1 \\ & \mathbf{A}\mathbf{x}\_t \le \mathbf{a} \qquad t = 0, 1, \dots, P\_H \\ & \mathbf{B}\boldsymbol{y}\_t \le \boldsymbol{\mathcal{B}} \qquad t = 0, 1, \dots, P\_H \\ & \mathbf{C}\boldsymbol{u}\_t \le \boldsymbol{\gamma} \qquad t = 0, 1, \dots, P\_H \\ & \mathbf{x}\_{sp} \in \mathbb{R}^{n\_{sp}} \end{aligned} \tag{4}$$

The notation adopted is the same as in the previous section, aside from the following: we treat the set-points together with the initial state of the system as uncertain parameters, thus resulting in a problem with simultaneous variations on the RHS and OFC. As indicated by (4), we consider both the initial states (*x*(*tk*)) and the various set-points (*xsp*) as uncertain

⎧

⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨

⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩

parameters. Notice that within an EWO framework, the deployment of such universal controllers is of grea<sup>t</sup> importance since they allow for rapid communication between the layer of control with RTO and scheduling. For instance, when integrating scheduling with control, the changeover times as well as the production rates are immediate results of the dynamic decisions made through the control system. In the next section, an algorithm for the design of such "multi set-point" explicit MPC is presented.

#### *3.2. Solution Algorithm for Analytical mp-NLPs with Global Uncertainty*

Here, we present an algorithm that can solve multi-parametric nonlinear programs with non transcendental nonlinear terms, i.e., nonlinear terms that have closed-form solutions. The proposed method can be seen as a generalisation of the our previous work on the solution of multi-parametric mixed integer polynomial programs [33] as well as the algorithm of Fotiou et al. [32] for multi-parametric polynomial programs. However, both of these methods were employed only for instances that the uncertainty is present on the right-hand side of the constraints and the latter does not compute the critical regions associated with each explicit solution.

The main idea of the algorithm proposed herein can be explained as follows: given a multi-parametric nonlinear program with analytical terms, or terms that can be expressed in a nontranscendental fashion, derive the first order KKT conditions and compute its solution using Gro¨bner bases by treating the uncertain parameters as symbols. The output of this step is a collection of candidate solutions which are explicit in *θ* and encompass: global and local optima as well as infeasible solutions. For these solutions, examine their dual and primal feasibility along with a constraint qualification so as to remove infeasible candidate solutions. Lastly, in order to report only the globally optimal solutions, perform a comparison procedure [33].

Problem (1) details a general formulation of multi-parametric programs. The case that *f* and/or **g** are analytically nonlinear and the uncertain parameters are in the OFCs along with the RHS and LHS of the constraints is used.

Deriving the 1st order KKT conditions of Problem (1) returns a system of equations that is square and is given by Equations (5) and (6).

$$\nabla\_{\mathbf{x}} L(\mathbf{x}, \boldsymbol{\theta}) = \mathbf{0} \tag{5}$$

$$
\lambda^T \mathbf{g} \left( \mathbf{x}, \theta \right) = \mathbf{0} \tag{6}
$$

*<sup>L</sup>*(**<sup>x</sup>**, *θ*) is the Lagrangian function of Problem (1), ∇**x** is the nabla operator with respect to the decision variables and *λ* are the Lagrange multipliers corresponding to the constraints. Because of the assumption that the nonlinearities have an analytical solution, Gro¨bner bases can be employed for the solution of the square system of equations because of its elimination property. Even though a tailored implementation of one of the already existing algorithms for computing Gro¨bner bases may be advantageous from a computational standpoint, it is beyond the scope of the present work and thus Mathematica 10 was employed as the computer algebra software in which the calculations were performed.

By solving the system of Equations (5) and (6), a number of candidate solution sets are returned. Note that although the original optimisation problem involves *nx* variables, in the current step, the variables for which we compute the explicit solution are *nx* + *ng*. The candidate solutions include the Lagrange multipliers together with the optimisation variables as explicit functions of the uncertain parameters, i.e., *λ*(*θ*) and **<sup>x</sup>**(*θ*), respectively.

#### **Definition 7.** *Candidate solutions [49]*

*A solution of the Problem (1) is said to be candidate if it satisfies the system of Equations (5) and (6) along with a constraint qualification, e.g., linear independence constraint qualification [15].*

In this part of the algorithm, due to strict complementary slackness, the collection of candidate solutions indicates the active and inactive constraints for each solution. Until this step, the set of solutions computed may be infeasible, local or global optima. By evaluating the primal and dual feasibility of the candidate solutions, the infeasible solutions can be rejected, i.e., Equations (7) and (8).

$$\mathbf{g}(\theta) \le \mathbf{0} \implies f \text{asibility conditions} \tag{7}$$

$$\lambda(\theta) \succeq \mathbf{0} \implies \text{optimality conditions} \tag{8}$$

Conditions (7) and (8) are evaluated by substituting the explicit expressions of the optimisation variables and they form a collection of parametric constraints. If for a candidate solution there exists a subset of the initial parametric space such that the aforementioned inequalities are satisfied, then this region is called the CR of the candidate feasible solution; otherwise, the candidate solution is infeasible and thus removed from further consideration. Note that the evaluation performed at this step, from a computer algebra perspective, is equivalent to computation of the corresponding CAD.

#### **Definition 8.** *Critical region [39,49]*

*A critical region (CR) is a partition of the parametric space where Conditions (7) and (8) are satisfied for a specific candidate solution. A critical region is characterised by a set of inactive/active constraints and can be discontinuous or nonconvex.*

In Algorithm 1, the pseudo-code of the presented method is given. The comparison procedure is outlined in [33].

#### *3.3. Illustrative Example*

The proposed methodology will be motivated through the following modified example by Domínguez et al. [29].

$$\min\_{x\_1, x\_2} x + 2x\_1^2 - 5x\_1 + x\_2^2 - 3\theta\_1 x\_2 - 6\theta\_2$$

Subject to:

$$\begin{aligned} 2x\_1 + x\_2 &\le 2.4 - \theta\_2 \\ 0.5\theta\_3 x\_1 + x\_2 &\le 1.5 \\ x\_1 \ge 0, \ x\_2 &\ge 0 \end{aligned} \tag{9}$$
 
$$\begin{aligned} \theta\_1 \le \theta, \ 0 \le \theta\_2 \le 4, \ 0 \le \theta\_3 \le 2 \end{aligned}$$

In the beginning, the first order KKT conditions of (9) are formulated and we derive the square system of Equations (10) and (14)

$$\nabla\_{\mathbf{x}} L(\mathbf{x}, \lambda, \theta) = \mathbf{0} \tag{10}$$

$$
\lambda\_1(2\mathbf{x}\_1 + \mathbf{x}\_2 - 2.5 + \theta\_2) = 0\tag{11}
$$

$$
\lambda\_2(0.5\theta\_3\mathbf{x}\_1 + \mathbf{x}\_2 - \mathbf{1}.5) = 0\tag{12}
$$

$$
\lambda\_3(-x\_1) = 0\tag{13}
$$

$$
\lambda\_4(-x\_2) = 0\tag{14}
$$

where *<sup>L</sup>*(**<sup>x</sup>**, *λ*, *θ*) is the Langangian of Problem (9). Equations (10) and (14) can be analytically solved through symbolic computations which return the dual and primal variables as functions of the uncertain parameters, i.e., *λ***(***θ***)** and *<sup>x</sup>***(***θ***)**. Systems (10)–(14) are solved in 0.006 s and, as shown in Table A1, fifteen candidate solutions are computed.

The primal and dual feasibility of the candidate solutions is examined by computing the CAD of the related disjunctions. If the result of this step is "False" the solution violates feasibility (in either the primal or dual sense), otherwise, a collection of explicit inequalities is returned which characterises the candidate solution's CR. The CAD computations of this example take 5.33 s and nine of them are nonempty. Despite the fact that nine candidate

0 ≤ solutions are primal and dual feasible, their global optimality is not guaranteed given the nonconvex nature of the problem. To this effect, for those regions, a new set of CAD computations is performed so as to identify overlaps.


It was found that *CR*10 and *CR*11 were overlapping, as shown by Figure 5, where the overlap (*CRint*) is shown as the dark partition in between the two CRs.

For the elimination of the resulting overlap, the comparison procedure is invoked and the logic disjunction, as illustrated below, is used for the CAD, as shown by Equations (15) and (16).

$$\exists \theta \mid \{\theta\_1, \theta\_2, \theta\_3\} \in \mathbb{C}R\_{int} \land z\_{\mathbb{C}R10}(\theta) \le z\_{\mathbb{C}R11}(\theta) \tag{15}$$

or

$$\left| \exists \theta \right| \left\{ \theta\_1, \theta\_2, \theta\_3 \right\} \in \mathcal{CR}\_{int} \land z\_{CR10}(\theta) \ge z\_{CR11}(\theta) \tag{16}$$

where *zCRi*(*θ*) denotes the optimal explicit value within *CRi*. The result of this step is partitioning of the parametric space where each explicit solution is the globally optimal. In this case, *CR*10 was shown to be dominant within the common parametric space and thus the overlap was subtracted from *CR*11. The algorithm terminates once no more overlapping critical regions are identified. In Table 2, an overview of the explicit solutions is presented, while in Figure 6, the final critical regions are shown. Practically, one would consult the CR column of Table 2 to identify based on the uncertain parameter values where the

uncertainty is realised, i.e., CR1 or CR2, and then the optimal cost can be computed by evaluating the corresponding expression from the "Explicit solution" column.

**Figure 5.** Instance of overlap between CRs.

**Figure 6.** Visualition of the critical regions for motivating example.


**Table 2.** Optimal explicit solutions of motivating example.

#### **4. Case Studies**

Here, we examine two case studies from process systems and their corresponding explicit controllers are designed. The computatonal experiments were performed on a workstation with 24 GB RAM, a 3.80 GHz processor and a Windows 10 64-bit operating system. For the symbolic calculation, the computer algebra system that was employed was Mathematica 10.2.

#### *4.1. Multiple-Input Multiple-Output Non-Isothermal CSTR*

We examine a non-isothermal MIMO multi-product CSTR where the decomposition reaction A →R happens under the kinetic law: −ℛ*<sup>b</sup>* = *krCb*. Additional details on the design and kinetics can be found in Camacho and Alba [50] and the data used for this case study can be found in Table 3. The system has two control inputs: the liquid (*Fl*) and coolant (*Fc*) flow rates, whereas the system's states are the temperature of the liquid (*Tl*) and the concentration of the decomposition product (*Cb*). Using the mass and energy balances, the dynamic model of the system is derived as given by Equations (17) and (18).

$$\frac{d(V\_l \mathbb{C}\_b)}{dt} = V\_l k\_l (\mathbb{C}\_{a0} - \mathbb{C}\_b) - F\_l \mathbb{C}\_b \tag{17}$$

$$\frac{d(V\_{l}\rho\_{l}\mathbb{C}\_{pl}T\_{l})}{dt} = F\_{l}\rho\_{l}\mathbb{C}\_{pl}T\_{l0} - F\_{l}\rho\_{l}\mathbb{C}\_{pl}T\_{l} + F\_{c}\rho\_{c}\mathbb{C}\_{pc}(T\_{c0} - T\_{c}) + V\_{l}k\_{r}(\mathbb{C}\_{a0} - \mathbb{C}\_{b})H \tag{18}$$

**Table 3.** Data of the multiple-input multiple-output CSTR case study.


Firstly, the systems (17) and (18) are transformed into an algebraic one in order to design the explicit controllers. Using the forward Euler method, the MPC problem of the discretised system is shown by Equations (19) and (20).

$$\min\_{\mathbf{u}} J(\boldsymbol{\theta}) = \sum\_{t=0}^{P\_H} \left\| \mathbf{x}(t) - \mathbf{x}\_{ref} \right\|\_2 \tag{19}$$

Subject to:

$$\begin{cases} \mathbb{C}\_{b\_{l+1}} = \mathbb{C}\_{b\_{l}} + \frac{h\_{c}(\mathbb{E}\_{l}(\mathbb{C}\_{a0} - \mathbb{C}\_{b\_{l}}) - \overline{\mathbb{V}}\_{l}\mathbb{C}\_{b\_{l}})}{V\_{l}}, \; 0 \le t \le P\_{H} - 1 \\ \mathbb{T}\_{l\_{l+1}} = T\_{l\_{l}} + h\_{c} \frac{\overline{F}\_{l\_{l}}\rho\_{l}\mathbb{C}\_{pl} - \overline{F}\_{l\_{l}}\rho\_{l}\mathbb{C}\_{pl}T\_{l\_{l}} + \overline{\mathbb{V}}\_{l}\rho\_{l}\mathbb{C}\_{pl}e\_{(T\_{l0} - \mathbb{C}\_{l})} + V\_{l}\mathbb{V}\_{l}(\mathbb{C}\_{a0} - \mathbb{C}\_{b\_{l}})}{V\_{l}p\_{l}\mathbb{C}\_{pl}} \\ 0.8 \le \mathbb{C}\_{b\_{l}} \le 3.5, \; 0 \le t \le P\_{H} \\ 280 \le T\_{l\_{l}} \le 400, \; 0 \le t \le P\_{H} \\ 0 \le F\_{l\_{l}} \le 1000, \; 0 \le t \le P\_{H} \\ 0 \le F\_{l\_{l}} \le 2000, \; 0 \le t \le P\_{H} \\ \mathbb{C}\_{b\_{l\_{l}:0}=0} = \theta\_{1}, T\_{l\_{l}:0=0} = \theta\_{2} \\ C\_{b\_{l\_{l}:0}=0} = \theta\_{3}, T\_{l\_{l}}^{ref} = \theta\_{4} \end{cases} (20)$$

By employing the presented solution algorithm for mp-NLPs, the related KKT system is solved analytically using Gro¨bner bases. It takes 0.76 s to compute 29 candidate solutions explicit in *θ*1, *θ*2 and a collection is shown in Table 4.



As mentioned in the previous section, the proposed algorithm can facilitate both continuous as well as discrete set-points for the solution of the resulting problem. For the MIMO case study, we consider 8 different set-points for which the explicit control law is derived. In Table A2, we outline the explicit solutions for the different set-points while the optimal partition of the uncertainty space is shown in Figure 7.

**Figure 7.** Critical regions for the mp-MPC controller for the MIMO CSTR.

We validate the performance of the explicit multi set-point controller by examining the transition between two steady states. We benchmark the controller's predictions against the globally optimal solution as computed using the BARON 14.4 solver. In Figure 8, the control and state evolution can be seen.

**Figure 8.** Comparative graphs of the control performance of the mp-MPC vs. NLMPC for case study 1.

#### *4.2. Isothermal Polymerisation CSTR*

Next, we examine the design of a multi set-point controller for the grade transition problem with polymerisation CSTRs. The model nonlinearities involve bilinear terms and square roots. The free radical polymerisation reaction happens in a CSTR that operates isothermally at 335 K, where methyl methacrylate (MMA) is produced [51,52]. The mathematical model is given by Equations (21) and (25). The system has 4 state variables, i.e., the concentrations of the monomer (*Cm*) and the initiator (*Cl*), the dead chains' molar concentration (*D*0) and the dead chains' mass concentration (*Dl*). The control input is the flow rate of the initiator (*Fl*) and one output, i.e., the molecular weight of the polymer produced (*y*). In the multiple steady states, different polymeric grades can be produced corresponding to alternative molecular weights. We provide the notation for this system in Table 5 and model parameter values can be found in Table 6.

$$\frac{d\mathbb{C}\_m}{dt} = -(k\_p + k\_{fm})\sqrt{\frac{2f^\*k\_l\mathbb{C}\_l}{k\_{Td} + k\_{Tc}}}\mathbb{C}\_m + \frac{F(\mathbb{C}\_{m\_{iw}} - \mathbb{C}\_m)}{V} \tag{21}$$

$$\frac{d\mathbf{C}\_{l}}{dt} = \frac{F\_{l}\mathbf{C}\_{l\_{in}} - F\mathbf{C}\_{l}}{V} - k\_{l}\mathbf{C}\_{l} \tag{22}$$

$$\frac{dD\_0}{dt} = (0.5k\_{T\varepsilon} + k\_{Td})\frac{2f^\*k\_l\mathcal{C}\_l}{k\_{Td} + k\_{T\varepsilon}}\mathcal{C}m + k\_{fm}\sqrt{\frac{2f^\*k\_l\mathcal{C}\_l}{k\_{Td} + k\_{T\varepsilon}}}\mathcal{C}\_{\text{ff}} - \frac{FD\_0}{V} \tag{23}$$

$$\frac{dD\_l}{dt} = M\_m(k\_p + k\_{fm})\sqrt{\frac{2f^\*k\_lC\_l}{k\_{Td} + k\_{Tc}}}C\_m - \frac{FD\_l}{V} \tag{24}$$

$$y = \frac{D\_l}{D\_0} \tag{25}$$

**Table 5.** MMA CSTR notation.



**Table 6.** Model parameters for the MMA polymerisation reactor.

The model nonlinearities are not transcedental and, thus, the presented method for the design of the multi set-point mp-MPC can be used. This polymerisation system has been examined intensively by the research community and it has been noted that online computation of its optimal control law can be challenging due to numerical instabilities arising because of scaling issues [52]. As a trade-off between computational complexity and stability of the integration scheme, we employ the forward Euler method with a step size of *h* = 36 *s*.

Following the proposed method, the globally optimal solutions are computed, whereas employing off-the-self global optimisation solvers for online implementation leads to extensive computational times. In Table 7, the results using the BARON 14.4 solver in GAMS with different prediction horizons are given.


**Table 7.** Computational effort for varying prediction horizons using BARON 14.4 [53].

† Reached time limit.

In this case study, eight uncertain parameters were considered, one for each state and set-point. In Table 8, we provide a mapping of the uncertain parameters of the control problem. Whilst having the set-point as continuous, uncertain parameters increase the computational complexity of the mp-P and one could argue that in the context of systems integration where the APC receives data by the real-time optimisation functionality, the same set-points may not always be realised. In such cases, following conventional mp-MPC frameworks, the explicit laws would have to be recomputed from the beginning (mp-P solution and implementation of the explicit solutions, possibly in a microchip), whereas, following the proposed framework, if the bounds of the set-points remain within the prespecified ranges, then the same multilayer controller can be readily used.


**Table 8.** MMA CSTR uncertain parameters.

The resulting mp-NLP involves one optimisation variable, eight uncertain parameters and ten constraints when a prediction horizon of unity is considered. Overall, we seek analytical solutions to eleven variables, i.e., the Lagrange multipliers and the optimisation variables. The solution of the mp-NLP returns five candidate solutions, as shown in Table A3.

Substituting the explicit expressions into the constraints, the feasible region is projected into the uncertainty space. Candidate solutions that satisfy primal/dual feasibility are considered for the next step of the algorithm; otherwise, they are discarded as infeasible. For instance, the 6th candidate solution violates dual feasibility as any value of *θ*6 would result in negative *λ*8.

Subsequently, the explicit inequalities for the remaining solutions are examined. The intersection of the feasible regions defined by the parametric inequalities defines the critical regions of the candidate solution. Because of the nonconvex nature of the problem, it is likely that explicit solutions may be valid in the same uncertainty space, thus overlapping. In order to compute only the global explicit solutions, we employ the comparison procedure. Three overlapping solutions were identified. An example of the inequalities defining the overlap between *CR*1, *CR*2 is shown by Equation (26).

$$\mathsf{CR}\_{ilt} \coloneqq \mathsf{CR}\_{1} \cap \mathsf{CR}\_{2} = \begin{cases} 0 \le \theta\_{3} \le 0.05 \\ 0 \le \theta\_{4} \le 300 \\ 0 \le \theta\_{5} \le 5 \\ 0 \le \theta \tau \le 0.05 \\ 0 \le \theta \tau \le 0.05 \\ 0 \le \theta\_{4} \le 300 \\ \theta\_{2} + 0.00242674 = 1.01114\theta\_{6} \\ 0 \le \theta\_{2} \le 0.199802 \end{cases} \qquad\qquad \begin{cases} 4.9899 \le \theta\_{1} \le 5 \\ 1617.74 + \frac{402802}{\theta\_{1}^{2}} \le \frac{16144.7}{\theta\_{1}} + \theta\_{2} \\ \theta\_{2} \le \frac{14397}{\theta\_{1}^{2}} \end{cases} \tag{26}$$

For illustration purposes, the mathematical definition of *CR*2 is given by Equation (A1) in the Appendix A. Due to the extensive set of inequalities defining the rest of the CRs, we do not detail them in the manuscript for the sake of space.

After the algorithm's convergence and with the optimal explicit solutions reported, the performance of the multi set-point mp-MPC's explicit control law is compared to that of conventional MPC. The solution of the online MPC is found by implementing the related NLP in GAMS and solving it to global optimality using BARON 14.4. As can be seen in Figure 9, the state and control evolution of the system are in perfect agreemen<sup>t</sup> when the two schemes are compared, thus highlighting the accuracy and correctness of the proposed framework while in Figure 10 the stability of the resulting control policy can be envisaged.

#### (**b**) *Fl* = *f*(*t*)

**Figure 9.** Plots comparing the control solutions computed by the proposed method (mp-MPC) and online NLMPC for the polymerisation CSTR.

**Figure 10.** Graph of the control cost function vs. time for the polymerisation CSTR.

We assess the performance of the controller for set-point tracking between two set-points. At the beginning, we assume the CSTR to be operated at a steady state of *y* = 15,000 *kg kmol* and then controlled towards *y* = 45,000 *kg kmol* where it is regulated for nine hours to produce a specific polymer grade. Next, the controller steers the system to the next set-point (*y* = 19,250 *kg kmol*), in which steady state another polymer is produced. The performance of the set-point tracking can be seen in Figure 11.

Finally, with respect to the scalability of the proposed method, a number of systems and prediction horizon settings were examined and, as shown in Table 9, for the current state of the art in computer algebra software, only small- to medium-scale systems can be efficiently facilitated.

**Figure 11.** Set-point tracking performance of the polymerisation reactor case study.


### **5. Conclusions**

We have presented a computer algebra-based algorithm for the analytical solution of mp-NLPs that involve uncertain parameters on the RHS and OFC as well as the LHS of the constraints. In the first step, Gro¨bner bases are used for symbolically expressing the optimisation variables and the Lagrange multipliers as functions of the uncertain parameters. Next, by computing cylindrical algebraic decompositions, the globally optimal CRs are defined. Building upon the proposed algorithm, we introduce a framework for the design of multi set-point explicit MPC for nonlinear systems. The proposed technique expands the scope of mp-MPCs, as we illustrate that it is feasible to design a single "multilayer" controller for capturing set-point tracking problems and potentially new model parameter estimations. Ongoing research focuses on the latter and how current progress in algebraic geometry can alleviate the related computational burden and allow for solutions of large-scale studies. Specifically, the application of machine learning techniques for faster evaluations of standard atomic formulas and thus reductions in the computational expense of CAD calculations is a promising direction.

**Author Contributions:** Conceptualization: V.M.C. and V.D.; Data curation: V.M.C.; Formal Analysis: V.M.C.; Funding acquisition: L.G.P. and V.D.; Investigation: V.M.C.; Methodology V.M.C., L.G.P. and V.D.; Visualization: V.M.C.; Writing—original draft: V.M.C.; Writing—review & editing: V.M.C. and V.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** The authors gratefully acknowledge financial support from EPSRC grants EP/M027856/1 and EP/M028240/1.

**Acknowledgments:** The authors gratefully acknowledge financial support from EPSRC grants EP/M027856/1 and EP/M028240/1.

**Conflicts of Interest:** The authors declare no conflict of interest.
