Multi Set-Point Explicit Model Predictive Control for Nonlinear Process Systems

In this article, we introduce a novel framework for the design of multi set-point nonlinear explicit controllers for process systems engineering problems where the set-points are treated as uncertain parameters simultaneously with the initial state of the dynamical system at each sampling instance. To this end, an algorithm for a special class of multi-parametric nonlinear programming problems with uncertain parameters on the right-hand side of the constraints and the cost coefficients of the objective function is presented. The algorithm is based on computed algebra methods for symbolic manipulation that enable an analytical solution of the optimality conditions of the underlying multi-parametric nonlinear program. A notable property of the presented algorithm is the computation of exact, in general nonconvex, critical regions that results in potentially great computational savings through a reduction in the number of convex approximate critical regions.


Introduction
High-fidelity and computationally efficient optimisation models are key for profitable decision making in process industries and have been the focus of extensive research over the years [1]. In recent years, the need for exploiting and explicitly considering interdependencies throughout the different layers of decision making has been underpinned by the enterprise-wide optimisation (EWO) concept [2]. Stemming from the progressively volatile and competitive market conditions, it is imperative for process industries to operate with agility in order to maximise their profitability [3]. EWO is aiming at increased profitability and resilience in process operations through the integration and simultaneous optimisation of existing information streams. Nonetheless, it comes at a considerable cost. Because of the multiple scales considered, EWO leads to computational challenges, thus preventing practitioners from harnessing the potential benefits such wide integration has to offer. Particularly, incorporating control considerations in an EWO fashion results in (mixed integer) nonconvex problems which are hard to solve.
By the same token, control considerations are ubiquitous in EWO problems. Figure 1 showcases how real-time optimisation and production scheduling exchange information with the layer of APC because of their interdependent decisions.
Real-time optimisation is concerned with the manipulation of systems' dynamics in order to achieve optimised profitability and operations. On the other hand, production scheduling determines the optimal allocation of resources for the completion of competing tasks. As indicated by Figure 1, both RTO and process scheduling exchange information with the layer of APC so as to achieve optimal dynamic operations. To this end, the research community has proposed different methods for their integration.
A common shortfall when focusing on integrating RTO and APC is that two different models are employed for the optimisation of the same system. Typically, a locally linear model of the initial nonlinear dynamics is used at the APC because of the need for fast A common shortfall when focusing on integrating RTO and APC is t two dierent models are employed for the optimisation of the same syst Typically, a locally linear model of the initial nonlinear dynamics is used the APC because of the need for fast solution rates while RTO considers original nonlinear model. This leads subsequently to issues related to subopti trajectories and non-reachable states (Young, 2006). Darby et al. [5], through their literature review regarding the integration of RTO and MPC, suggested that for a successful integration, common issues such as model mismatch among the layers of RTO and APC should be eliminated. Nonetheless, in real industrial processes, model degradation and other factors can result in model mismatch, so the consideration of parameter estimation and data reconciliation functionalities is needed to integrate RTO and MPC, as indicated by Figure 2. The interaction between real-time optimisation and model predictive control can be categorised broadly into three classes: (i) dynamic RTO (d-RTO), (ii) static RTO (s-RTO) and (iii) economic model predictive control (e-MPC). Both s-RTO and d-RTO are twolayer schemes where reference trajectories are passed to the layer of APC in the form of set-points [6]. While under the static real-time optimisation paradigm, the optimisation problem is solved at specific instances whenever new data become available or when steady state is achieved, in the d-RTO paradigm, the system's transient behaviour is explicitly considered, thus resulting in dynamic optimisation problems. e-MPC [7] refers to singlelayer strategies which are incorporated into the control structure economic considerations. In that spirit, De Souza et al. [8] proposed the inclusion of the gradient of the economic objective function into the MPC objective as a single-layer strategy. Considering uncertain systems, Chachuat et al. [9] examined alternative model adaptation strategies. This article is motivated by the abovementioned issues and aims at introducing a method for designing multi set-point explicit controllers for nonlinear systems through recent advances in multi-parametric programming. Multi-parametric programming (mp-P) has received considerable attention from the process systems engineering community because of its unique ability to aid in the design of explicit model predictive controllers and thus shift the computational burden associated with offline control [10]. We examine a case of multi-parametric nonlinear programs (mp-NLPs) that involve both endogenous uncertainty, in the form of left-hand side parameters (LHS), as well as exogenous uncertainty in the cost coefficient of the objective function (OFC), and, on the right-hand side of the constraints (RHS), uncertain parameters on the right-hand side (RHS). In engineering problems, LHS uncertainty arises from variations in model coefficients, due to parameter estimation errors or model mismatch; OFC uncertainty arises due to fluctuation in market prices or control penalties while RHS uncertainty can be due to varying system exogenous factors. The contribution of the present work is a novel framework for the design of multi set-point explicit controllers for nonlinear process systems.
The remainder of the article is organised as follows: in Section 2, an overview of the field of multi-parametric programming and explicit MPC is given, and then, in Section 3, the proposed algorithm is detailed and a framework for multi set-point explicit controllers is introduced. In Section 4, two case studies are examined so as to illustrate the main computational steps of the proposed methodology and, lastly, in Section 5, concluding remarks are made.

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: where θ stands for the vector of uncertain parameters, which is n θ -dimensional, x is the n x -dimensional vector of continuous decision variables, g(x, θ) is the vector of inequality constraints and f is the objective function as a mapping R n x ×n θ → R, both of which are assumed to be C 2 (twice continuously differentiable). Problem (1) is a multi-parametric program and its solution results in the partition of the parametric R n θ -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 z(θ), 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 n θ 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.

mp-NLP Solution Techniques RHS LHS OFC Comments
Dua and Pistikopoulos [21] --Convex Johansen [23] --Convex Acevedo and Salgueiro [22] --Convex Johansen [25] --Convex Dua et al. [28] --Nonconvex Fotiou et al. [31] --Polynomial Narciso [26] --Convex Domínguez and Pistikopoulos [24] --Convex Charitopoulos and Dua [33] --Polynomial Pappas et al. [34] --QCQPs 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): where x t , u t , z t are the state, control input and system output vectors, respectively, at every sampling point, t, and are n x , n u , n y -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 n x +n u → R is a stage cost and E : R n x → R is a terminal cost function. The repetitive solution of Problem (2) provides the optimal cost Φ(x(t k )) and the optimisation vector, i.e., the sequence of optimal control inputs u * = u * 1 , u * 2 , . . . , u * P H −1 over the finite prediction horizon P H . 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).
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. 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.

Gröbner 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.
Grö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, Grö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 Gröbner basis is, it would be useful to define some preliminary concepts.

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

Definition 2. Term order
A term order is defined with regard to a set of power products (T n ) 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 Gröbner bases theory.

Definition 3. Ideals
Let R be a field and R[x 1 , x 2 , . . . , x n ] be a ring over the field of n-variate polynomials. Let a finite subset of the field, G = {g 1 , g 2 , . . . , g t }, then an ideal I can be generated by G as follows: For problems that accept analytical solutions, their existence is guaranteed by the Hilbert basis theorem [43], which also guarantees that algorithms used to compute Gröbner bases can terminate in a finite number of steps.

Definition 4. Gröbner basis [41]
A set of non-zero polynomials G = {g 1 , . . . , g t }, contained in an ideal I, is called a Gröbner basis for I if and only if for all g ∈ I, such that g = 0, there exists i ∈ {1, . . . , t} such that l p(g i ) divides l p(g), where l p(· ) stands for the leading power product of a polynomial function.
For the calculation of Gröbner bases, apart from Buchberger's algorithm, Faugè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 Gröbner bases to be created. Lastly, software implementations of different algorithms that compute Grö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.

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.
Let R[x 1 , x 2 , . . . , x n ] 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: 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: ). Semi-algebraic sets of R n can be written as a finite union of semi-algebraic sets of the form: where g 1 . . . , g ω , g ω+1 , . . . , g t are in g ∈ R[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].

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 great 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 n u set-point, is to solve n u mp-P problems and thus design n u 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. The design problem of a multi set-point mp-MPC can be formulated as in (4). The design problem of a multi set-point mp-MPC can be formulated as in (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(t k )) and the various set-points (x sp ) as uncertain parameters. Notice that within an EWO framework, the deployment of such universal controllers is of great 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.

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 Grö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).
L(x, θ) 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, Grö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 Grö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 n x variables, in the current step, the variables for which we compute the explicit solution are n x + n g . The candidate solutions include the Lagrange multipliers together with the optimisation variables as explicit functions of the uncertain parameters, i.e., λ(θ) and x(θ), 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).
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. [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.

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

Illustrative Example
The proposed methodology will be motivated through the following modified example by Domínguez et al. [29].
Subject to: In the beginning, the first order KKT conditions of (9) are formulated and we derive the square system of Equations (10) and (14) ∇ x L(x, λ, θ) = 0 (10) where L(x, λ, θ) 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 x(θ). 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 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.  Evaluate with a first order constraint qualification, e.g. Linear Independence Constraint Qualification (LICQ) 9 Evaluate with primal and dual feasibility conditions (Cylindrical Algebraic Decomposition computation): if CR = ∅ then Candidate solution i is infeasible and discard from TEMP. 12 else 13 Candidate solution i is feasible and append CR i to TEMP. It was found that CR 10 and CR 11 were overlapping, as shown by Figure 5, where the overlap (CR int ) 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).
where z CRi (θ) denotes the optimal explicit value within CR i . 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.
and nine of them are non-empty. Despite the fact that nine candidate solutions are primal and dual feasible, their global optimality is not guaranteed given the non-convex nature of the problem. To this eect, 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 Fig. 5 where the overlap (CR int ) is shown as the dark partition in between the two CRs.

CRs
Explicit Solution

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.

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: −ℛ b = k r C b . 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 (F l ) and coolant (F c ) flow rates, whereas the system's states are the temperature of the liquid (T l ) and the concentration of the decomposition product (C b ). Using the mass and energy balances, the dynamic model of the system is derived as given by Equations (17) and (18).  (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).
Subject to: By employing the presented solution algorithm for mp-NLPs, the related KKT system is solved analytically using Grö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. 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.

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 (C m ) and the initiator (C l ), the dead chains' molar concentration (D 0 ) and the dead chains' mass concentration (D l ). The control input is the flow rate of the initiator (F l ) 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.    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. 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.

Parameter
Range Notation 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). (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 agreement 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.  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.

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, Grö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.  Table A1. Candidate solutions of motivating example.