Design of multi-parametric NCO tracking controllers for linear dynamic systems

A methodology for combining multi-parametric programming and NCO tracking is presented in the case of linear dynamic systems. The resulting parametric controllers consist of (potentially nonlinear) feedback laws for tracking optimality conditions by exploiting the underlying optimal control switching structure. Compared to the classical multi-parametric MPC controller, this approach leads to a reduction in the number of critical regions. It calls for the solution of more di ﬃ cult parametric optimization problem with linear di ﬀ erential equations embedded, whose critical regions are potentially nonconvex. Examples of constrained linear quadratic optimal control problems with parametric uncertainty are presented to illustrate the approach.


Introduction
Driven by the need to improve performance and reduce economic costs in industrial processes, on-line optimization and real-time control have been receiving a lot of attention. Many such industrial applications involve fast dynamic systems operated under constraints, typically reflecting physical operation bounds and/or safety requirements [1]. Optimal control strategies can be determined by solving constrained optimization problems based on a dynamic model of the system. One major challenge with this approach is how to effectively deal with uncertainty stemming from model mismatch and process disturbances, as optimal operation needs to be decided on-line using real-time feedback. The strategy employed in classical model predictive control (MPC) [2] handles this uncertainty by repeatedly solving the optimization problem on-line in order to update the optimal inputs. This is often a rather computationally demanding task that may cause serious delays especially for systems with fast dynamics, leading to suboptimal performance or even infeasibility.
Several approaches have been proposed in the literature to overcome the need for repeatedly solving optimization problems on-line. In the multi-parametric programming paradigm [3], the optimization is performed off-line, resulting in a priori explicit mapping of the solutions, effectively control strategies, as a function of measurable quantities. For continuous-time systems, this approach calls for the solution of multi-parametric dynamic optimization problems (mp-DO) [4].
In practice, such mp-DO can either be transformed into a finite-dimensional multi-parametric program via control vector parameterization [4] or handled directly into its native infinite-dimensional form using the corresponding optimality conditions. In the context of solving the finite-dimensional multi-parametric program, there exist numerous publications [3,[5][6][7][8][9][10], whereas solving the infinite-dimensional counterpart has received relatively little attention [4].
Parametric optimal solutions for the latter can be obtained by sensitivity analysis, also known as neighboring extremal (NE) control. In this approach a feedback law is derived in the neighborhood of a nominal solution, where the switching structure-namely, the sequence of active path and terminal constraintsremains the same [11][12][13][14]; see also [12,15,16] for a discussion of the differentiability and stability of parametric solutions. For constrained linear-quadratic optimal control problems in particular, Sakizlis and coauthors [4] have shown that the mp-DO can be written as a multi-parametric boundary value problem using the Pontryagin Minimum Principle [17][18][19]. The continuous-time optimal control trajectory is expressed as a time-varying functions of selected parameters, which provides a means for determining the control switching structure using standard multi-parametric programming techniques.
Another approach to reducing the on-line computational burden involves tracking the necessary conditions for optimality, namely NCO-tracking [20,21]. There, an optimal control policy is obtained indirectly by forcing the NCOs to zero. Such a forcing requires knowledge of the switching structure of the optimal control, based on which feedback control laws can be derived on account of the available output measurements, effectively converting an optimal control problem into a set of self-optimizing feedback control laws [20,22]. However, a key assumption for this controller to enable optimal or even feasible operation is that the switching structure should remain unchanged, which may not be the case when large uncertainty is present.
This paper presents a methodology for combining mp-DO and NCO-tracking into a unified framework for model predictive control, for which we coin the name multi-parametric NCO-tracking control. Such a combination is especially promising in that multi-parametric programming provides a means for relaxing the fixed switching structure assumption in NCO-tracking, thereby paving the way towards a theoretical justification for NCO-tracking too. In addition, the use of feedback laws tracking the optimality conditions inside multi-parametric controller provides a means for reducing the number of critical regions compared with classical mp-MPC controllers based on control vector parameterization.
The rest of the paper is organized as follows. Sect. 2 provides some background for both multiparametric programming and NCO tracking method. Sect. 3 presents the multi-parametric NCO-tracking methodology, including the use of mp-DO for mapping subregions of the uncertain parameters to various switching structure and the implementation of the corresponding NCO-tracking controllers in a receding horizon manner. Several numerical examples are given in Sect. 4 to illustrate the approach. Finally, Sect. 5 concludes the paper.

Multi-parametric MPC
Multi-parametric programming provides a means for computing the solution mapping of an optimization-based control problem off-line, based on a model of the dynamic system. The optimal control trajectory is expressed as a function of given parameters, usually some uncertain measured quantities, thus avoiding the need for repeatedly solving optimization problems on-line when these parameters vary [5,6,9,23,24]. In the context of MPC, this approach is referred to as mp-MPC or explicit MPC. A comparison of the control framework of MPC and mp-MPC is shown in Fig. 1. Instead of repeatedly solving the optimization problem during run-time as in MPC, mp-MPC computes a mapping between the uncertain parameters and their corresponding optimal solutions off-line, then simply selects the pre-computed control law at run-time after the uncertainty is revealed. In the multi-parametric solution, the parameter space is partitioned into a number of critical regions, where the optimal input variable, u, is described by a given function of the parameters, θ, as Each critical region corresponds to a unique combination of active constraints, the boundary of which can be computed from sensitivity analysis of the KKT conditions by keeping the inactive constraints non-binding and the multipliers associated with the active constraints non-negative. A basic procedure for determining the critical regions is the following [6]: 0. Define the uncertainty domain Θ, and set N = 0. 1. Select a feasible point θ in the region Θ \ ∪ N i=1 CR (i) . If no such point exists, stop; else, set N ← N + 1. 2. Construct the critical region CR (N) around θ, wherein the active constraints are the same, e.g. using sensitivity analysis of the KKT conditions. 3. Return to step 1. 4. Unify the regions and solutions for a more compact representation.
A linear mp-MPC problem is considered next to illustrate this procedure. The same example will be used throughout the theoretical part of the paper to illustrate the developments.
Illustrative example. We consider the problem to steer the state x(t) of the following second-order system to zero, by manipulating the bounded input u(t) ∈ [−2, 2]: The mp-MPC problem obtained by discretizing the dynamics on N t time intervals along the time horizon where the parameters θ ∈ [−2, 2] 2 corresponds to the initial state; the matrices F x and F u in the discretized are given by with the sampling time T := 1/N t ; and the weighting used in the objective function is as follows Numerical solution of the mp-QP (3), here using the PAROC framework [8], provides expressions of the optimal controls u = [u(0), . . . , u(N t − 1)] as explicit functions of the initial conditions θ, in the form (1). The critical regions for the optimal solution are shown in Fig. 2(a) in the case N t = 10, whereby each region CR i corresponds to a piecewise affine functions u = K i θ + k i , with K i ∈ R N×2 and k i ∈ R 2 . Here, the region labeled CR08 in the central part corresponds to the case that none of the input constraints are active. The regions above CR08 correspond to the input lower bound being active for one or more time intervals, and the farther from CR08 the larger the number of active constraints. Likewise, the regions below CR08 correspond to the input upper bound being active for one or more time intervals.
For comparison, critical regions in the case N t = 20 are shown in Fig. 2(b). Here again, the region in the center corresponds to unconstrained solution and other regions, either above or below it, correspond to input constraints being active for the first few time intervals. The multi-parametric solution becomes more accurate due to the use of a smaller sampling time, but at the same time the number of critical regions increases significantly, thereby defining a trade-off between accuracy and computational tractability. In contrast, the approach proposed in this paper removes the need for discretizing the dynamics and the control trajectories, in order to reduce the number of critical regions.

NCO tracking
NCO tracking is a measurement-based optimization approach to enforcing optimality in the presence of uncertainty, via tracking of the necessary conditions for optimality (NCO). This way, a dynamic optimiza-  tion problem is transformed into a feedback control problem [20,25], which may lead to a large reduction of the on-line computational effort by avoiding the repeated solution of an optimal control problem.
The design of the NCO-tracking controllers starts by detecting the switching structure of the optimal control in order to formulate a feedback strategy via appropriate pairing between the input and output variables-the so-called solution model [26,27], see Fig. 3(a). Along each arc, a certain combination of inputs may be used for tracking the active path constraints, whereas the remaining inputs are adapted for forcing stationarity conditions (gradients) to zero. This latter forcing usually calls for approximation techniques, such as neighboring-extremal control [28][29][30][31] or extremum-seeking control [32,33]. It is sometimes possible to arrive at a fully decentralized control scheme, for instance using directional information [20,22,34].
A key limitation with the current NCO-tracking methodology nonetheless lies in the fact that the underlying optimal control switching structure might change in the presence of uncertainty. As a result, the NCO-tracking controller may be suboptimal and could even lead to infeasible operation due to constraint violation. Although the assumption of a constant structure is often satisfied in batch process optimization applications [35], it is not well suited for MPC applications where constraints frequently activate or deactivate. It has been suggested that the control switching structure could be monitored by some supervisory system [28]. The developments in Sect. 3 provide the foundations for such an approach to handling a varying optimal switching structure. It relies on the mapping between uncertain parameters and optimal switching structures using mp-DO, and the subsequent formulation of optimal control laws that can be applied in a receding horizon manner, namely multi-parametric NCO-tracking controllers.

Problem statement
The main contribution of the paper is a methodology for the derivation of multi-parametric NCOtracking controllers for constrained linear-quadratic optimal control problem in the form: where Φ is the optimal value function; x(t) ∈ R n x and u(t) ∈ R n u are the state variables and input variables, respectively, at a given time t; t 0 and t f are the initial and final times; g : R n x × R n u → R n g and h : R n x → R n h define the path and terminal inequality constraints; and Q f 0, Q 0, and R 0 are given weighting matrices. We assume that the uncertain parameters θ ∈ R n θ appear linearly in the initial conditions, dynamics, path constraints, and terminal constraints of problem (P θ ).
The proposed methodology involves two steps, as depicted in Fig. 3(b): • The first (off-line) step defines the multi-parametric control structure, namely mapping the optimal control structure to given measurable quantities, such as the uncertain initial conditions, using mp-DO. This results in a partitioning of the uncertain parameter domain into a number of critical regions, each corresponding to a unique sequence of active path constraints and active terminal constraints. As well as giving a set of conditions for characterizing each critical region, mp-DO determines feedback control laws in the form: where the junction times t (i) 1 (·), . . . , t (i) (·) in the optimal switching structure for critical region CR i are themselves dependent on θ.
• In the subsequent (on-line) step, the multi-parametric NCO-tracking controller is applied in a receding horizon manner. This involves determining the critical regions containing the measured parameters θ and applying the corresponding feedback law until a new measurement becomes available at the next sampling time. Because the switching time functions t (i) k (·) are typically defined implicitly in practice, even for constrained linear-quadratic control problems, one can either derive fully explicit feedback laws by approximating this functional dependency, or else apply a Newton iteration to compute the t (i) k for given values of θ at each sampling time. Both steps are detailed subsequently.
Notation. D i x f denotes the i-th partial derivative of a function f with respect to x and f ( j) , the j-th order derivative of with respect to t. The path constraint g i is said to be of order (or degree) For simplicity, we introduce the notation and g 0 , for each j = 1 . . . σ i . We also make use of the notation Finally, by a slight abuse of the notation, an overbar is used to indicate subsets of the terminal or path constraints that are active along a given arc, such asḡ(x(t), u(t)) =Ḡ x x(t) +Ḡ u u(t) +Ḡ θ θ +ḡ 0 ≤ 0 andμ(t).

Multi-parametric dynamic optimization 3.2.1. Solution structure
For each instance of the parameters θ, the optimal solution of problem (P θ ) exhibits a certain switching structure, denoted by S(θ). The sequence of active path constraints and active terminal constraints can be characterized by solving the first-order NCO for Problem (P θ ), which come in the form of a multipoint boundary value problem [17]. Under the assumption that the number of arcs N t (θ) is finite for each parameter value [36], we denote by t k (θ), k = 1 . . . N t − 1(θ), the sequence of junction times for each arc in S(θ), with t 0 (θ) = t 0 and t N t (θ) = t f . These times correspond to the activation or deactivation a particular path constraint or to a touch-and-go point for a higher order path constraint. We denote the sets of path constraints activating, deactivating, or contacting at t k (θ) by EN k (θ), EX k (θ) and CO k (θ), respectively. Moreover, AC k (θ) and NAC k (θ) denote the sets of active/inactive path constraints along the kth arc, ; and AC f (θ) and NAC f (θ), the sets of active/inactive terminal constraints. 7 Besides its switching structure, characterizing an optimal solution involves determining: (i) the quadruplet of trajectories (u(t), x(t), p(t), µ(t)) along each arc, where p(t) ∈ R n x are the co-state (adjoint) variables, and µ(t) ∈ R n g , the multipliers for the path constraints; (ii) the values of the multipliers ν ∈ R n h for the terminal constraints; and (iii) the values for the multipliers π j k,i for j = 1 . . . σ i − 1, i = 1 . . . n g , k = 1 . . . N t (θ) − 1 at points of discontinuity of the co-state trajectories p(t) (if any). Provided certain controllability and regularity assumptions hold (see below), the following conditions must be satisfied at an optimal solution of Problem (P θ ), according to the indirect adjoining approach [19]: for each i = 1 . . . n g and each j = 1 . . . σ i , and with the Hamiltonian function : for each i = 1 . . . n g and each j = 1 . . . σ i .
Note that the multipliers π may only appear in the optimal conditions for problems with pure state path constraints of order 1 or higher; they can be discarded in problems having mixed control-state path constraints only, where the adjoint trajectories are continuous at junction times. In general, the foregoing optimality conditions (5)- (17) are only necessary under the additional assumptions that: (i) the pair (F x , F u ) is controllable, which precludes abnormality [37]; and (ii) both the active path and active terminal constraints are regular [19], Moreover, under the extra assumption of strict complementarity slackness for the multipliers ν, π j k,i and µ(t) for a given parameter value θ 0 , and by strict convexity of the objective function and linearity of the dynamics and constraints, the optimal trajectories u(t), , optimal multipliers ν and π j k,i , and optimal switching/contact times t k describe differentiable functions in an (open) neighborhood of θ 0 [14][15][16]; see also [4]. Expressions for these functions are established in the following subsection.

Feedback control laws
Using the previous optimality conditions, explicit feedback control laws can be derived along each arc of the optimal solution. Using condition (7), together with the fact thatḡ (σ) θ θ +ḡ (σ) 0 = 0 along an arc, we havē which are both well-defined under the assumption thatḠ (σ) u is full rank. In turn, the state and co-state equations (5)-(6) can be rewritten in the form with: This way, we may express x(t) and p(t) at each t ∈ [t + k−1 (θ), t k (θ) − ], and therefore also u(t) and µ(t), as parametric functions of the uncertainty θ, the junction times t k , and the state/co-state values at t k : In the case that A (k) xp is nonsingular, we have: When A (k) xp is singular, an explicit expression can be obtained by considering the normal Jordan form of A (k) xp instead. At this point, parametric expressions for the terminal and interior-point constraint multipliers ν and π k can be obtained by piecing together (21) on [t 0 , t f ] and exploiting the equality conditions in (10)-(11) and (14)- (17). In the case of mixed state-input path constraints only, the optimal state and co-state trajectories are both continuous at the junction times. Then, sinceH x is full rank, and provided that A (k) xp is invertible on each arc, parametric expressions for the active terminal constraint multipliersν, terminal state x(t f ) and initial adjoint p(t 0 ) can be obtained from the following linear system: (In the case of a single control arc, N t (θ) = 1, the term N t (θ) j=k+1 A ( j) xp [t j (θ) − t j−1 (θ)] in the right-hand side of (24) cancels to zero.) Overall, for a given structure S(θ), the solution of Problem (P θ ) can therefore be expressed in parametric form as Naturally, this construction can be automated in a practical implementation. One could also exploit the remaining optimality conditions (13) in order to determine parametric expressions of the junction times t k as a function of θ alone. Explicit expressions are often not possible, however, due to the inherent nonlinearity of the parametric state/co-state expressions (21) in t k and θ. In practice, one may either use approximate explicit expressions for t k (θ), or compute these junction times on-line using a Newton iteration. These considerations are discussed further in Sect. 3.3.

Critical regions
Each critical region corresponds to a subset of the uncertain parameter domain Θ ⊆ R n θ , whereby the corresponding optimal control solutions all share the same switching structure in terms of the sequence of active path constraints and active terminal constraints. Formally, given the switching structure S comprised of N t arcs with corresponding index sets EN k , EX k , CO k , AC k , NAC k , AC f , and NAC f , the critical region CR S associated with S is defined as The boundary between two critical regions thus corresponds to either an inactive terminal or path constraint activating, or an active terminal or path constraint inactivating, or two subsequent junction times becoming equal. Similar to the procedure outlined in Sect. 2.1 for the construction of critical regions in mp-MPC, a systematic procedure for constructing critical regions in mp-DO is as follows: 0. Define the uncertainty domain Θ ⊂ R n θ , and set N = 0. 1. Select a point θ ∈ Θ \ ∪ N i=1 CR (i) that is feasible for (P θ ). If no such point exists, stop; else, set N ← N + 1. On termination, this procedure returns the number N of critical regions contained in the initial domain Θ, together with the solution structure S (i) and the feedback laws F (i) in each region CR (i) , i = 1 . . . N, as defined implicitly by (26). A number of remarks are in order regarding the practical implementation of this procedure: • The selection of a point θ ∈ Θ \ ∪ N i=1 CR (i) in step 1 is a difficult problem in general, since (26) may describe non-convex subsets. Suppose, without loss of generality, that the critical regions can be written in the form CR In particular, techniques based on complete search have become available in recent years to address such constraint satisfaction problems; see, e.g., [38]. More and more constraints are appended as the number of critical regions found increases, and the problem eventually becomes infeasible when the parameter region has been exhausted.
• Step 2 involves solving the constrained dynamic optimization problem (P θ ) and characterizing its underlying solution structure. A numerical solution can be computed using direct solution methods [39], which discretize the control and/or state variables in order to arrive at an approximate, finitedimensional optimization problem. In the present case, this approximate problem is a convex QP for which efficient and reliable large-scale solvers are available. Moreover, such direct solution methods can be used in combination with adaptive structure detection techniques, as discussed e.g. in [26,40].
The computation of critical regions and characterization of the corresponding feedback laws is now illustrated for the linear mp-MPC problem first considered in Sect. 2.1.
Illustrative example (continued). We consider a mp-DO problem for the second-order dynamic system (2) in the form of Problem (P θ ), with and F θ , f 0 , G x , G θ , H x , H θ and h 0 all zero. This mp-DO aims to solve the same optimization-based control problem as in (3), yet in continuous-time form-that is, without discretizing the time horizon and the dynamics. The solution for Θ = [−2, 2] 2 using the approach described in this subsection yields N = 3 critical regions, as shown in Fig. 4(a): • The optimal control strategy for an initial state in the critical region CR01 is comprised of a unique unconstrained arc, which is therefore identical to classical unconstrained LQR; • For an initial state in CR02, the optimal solution structure is comprised of two arcs, namely a constrained arc with the input at its lower bound followed by an unconstrained arc; • For an initial state in CR03 likewise, the optimal solution consists of a constrained arc with the input at its upper bound followed by an unconstrained arc.
Expressions for the feedback laws F (i) in each region CR (i) , i = 1 . . . 3 are reported in Appendix A. For instance, the open-loop optimal trajectories for the initial states θ (1) = (1, −0.8) ∈ CR01 and θ (2) = (1, 1) ∈ CR02 are shown in Fig. 4(b) and 4(c), respectively. Note that an explicit expression for the switching time t 1 in both CR02 and CR03 as a function of the initial conditions θ can also be found for this simple problem The critical regions of the mp-DO problem are closely related to those of the discretized mp-MPC problem in Fig. 2. First of all, the central region CR01 in Fig. 4(a) matches the central regions CR08 and CR15 in Fig. 2(a) and 2(b), respectively, and all three feedback laws correspond to the same unconstrained LQR strategy, u(t) = [8.198 8.198]x(t). The rest of the critical regions in Fig. 2 have their first few control stages either at the upper bound or at the lower bound, and the farther the central region, the larger the number of time intervals with active constraints. This behavior is fully consistent with the switching time t 1 between the upper/lower bound and the interior arc in the continuous-time optimal control increasing in (a) Critical regions moving further away from CR01. By construction, the feedback laws corresponding to the critical regions other than the central region in Fig. 2 are approximations of the feedback laws computed for the critical regions CR02 and CR03 in Fig. 4(a). The finer the discretization, the closer the approximate feedback laws to the continuous-time ones; but this comes at the price of a significant increase in the number of critical regions in order to better capture the variations in switching time, about double the number from N = 10 to N = 20.
In this interpretation, mp-DO thus provides a means of reducing the number of critical regions by accounting for the underlying switching structure and introducing the switching times in the feedback law parameterization. This reduction becomes more effective as certain arcs in the mp-DO solution spans many time intervals in the discretized mp-MPC problem.

Multi-parametric NCO-tracking controller
The critical regions CR (1) , . . . , CR (N) and their corresponding feedback laws F (1) , . . . , F (N) define a multi-parametric NCO-tracking controller for the problem (P θ ) at hand. At this point, we like to reiterate that the feedback laws are determined off-line via a (continuous-time) mp-DO formation, whereas the closed-loop mp-NCO-tracking controller selects the appropriate control law after the uncertainty has been revealed and applies them in a receding horizon manner. In this respect, mp-QP and mp-MPC are the discrete-time counterparts of mp-DO and mp-NCO-tracking, respectively.
The application of the mp-NCO-tracking controller, in a receding horizon manner, involves determining the critical region corresponding to the uncertainty revealed at the current sampling time. Checking whether or not θ ∈ CR (i) for a given i = 1 . . . N can be done in two steps: 1. Given the feedback laws F (i) , determine the switching times t 1 (θ), . . . , t N (i) t (θ)−1 (θ) satisfying the Hamiltonian continuity conditions (13); 2. Test whether all the auxiliary inequalities defining CR (i) as in (26) are satisfied.
Once the correct critical region CR * has been identified, the controller simply applies the associated feedback law F * until the following sampling time.
In this approach, the switching times may be computed using a Newton iteration (or a robustified variant such as the dog-leg method). The on-line computational burden can be reduced by initializing the switching times with the values at the previous sampling time. This warm-starting strategy is most effective when the sampling frequency is fast, so the variations in switching times between two sampling times remain small. Moreover, since only the first part of the optimal feedback law is applied in a receding horizon strategy, the same control law may be applied several times when not foreseeing any switching time nearby. Detection of such instances can lead to large computational-time reductions. An alternative approach involves approximating the functional dependency of the switching times with respect to the uncertainty, for instance using linear or polynomial functions in θ. These extra laws can be determined via parametric programming too, possibly after further partitioning of the critical regions for keeping the approximation level under control.
The application of a multi-parametric NCO-tracking controller for the same linear mp-MPC problem as in Sect. 2.1 is presented below.
Illustrative example (continued). The critical regions and feedback laws given in Appendix A and shown in Fig. 4(a) define a multi-parametric NCO-tracking controller for the linear system (2). The application of this controller in a receding horizon scheme is illustrated in Fig. 5, for a sampling period of ∆T = 10 −2 and from the initial state θ (2) = (1, 1) ∈ CR02. The control and state trajectories in the noise-free case ( Fig. 5(a)) closely match those of the nominal, open-loop case in Fig. 4(c). A difference here is the value of the switching time t 1 , which is postponed with the closed-loop controller due to the finite sampling period. The closed-loop trajectories for the same problem, now with Gaussian white noise with signal-to-noise ratio of 50dB added to the state measurements at each sample time, are shown in Fig. 5(b) for comparison. Recall also that explicit expressions of the switching time t 1 as a function of θ can be obtained in both critical regions CR02 and CR03, so no Newton iteration or approximation is needed in this instance.

Computational aspects
The computational burden for the proposed multi-parametric NCO-tracking framework involves two distinct components, namely the off-line controller design and its on-line application. The former is dominated by the solution of an mp-DO problem, where all the variables are linear functions of the parameters, but the switching times. Because of this nonlinearity, the critical regions do not describe convex polytopes anymore, but yield general non-convex regions. Even though the explicit characterization of a region's boundary is not needed, the enumeration procedure in Sect. 3.2.3 calls for the application of complete search methods in order to find a new critical region or to establish that all the critical regions have been mapped already, which can be computationally demanding (if at all tractable) [38]. In practice, one can apply model reduction techniques for reducing the order of the dynamic model and improve the computational tractability of the mp-DO problem, yet without causing too big a performance loss for the resulting controllers.
With regards to the on-line application aspects, mp-NCO-tracking controllers provably reduce the number of critical region compared with their discretized mp-QP counterparts. Nonetheless, the on-line computational burden with mp-NCO-tracking is typically dominated by the determination of the switching times corresponding to the measured uncertainty θ, e.g. using a Newton iteration. In practice, efficient warmstarting strategies could be developed for minimizing this burden, especially for short sampling periods. Whether or not such strategies will make mp-NCO-tracking competitive with simple look-up table evaluation in discretized mp-QP despite a possibly much larger number of critical regions will be explored as part of future work.

Numerical case studies
The objective of the numerical case studies in this section is two-fold: (i) illustrate the computation of mp-DO critical regions for problems with first-and higher-order state path constraints along with terminal constraints, thereby complementing the foregoing illustrative example with simple bound constraints only (Sects. 4.1 and 4.2); and (ii) present the mp-DO construction and the corresponding multi-parametric NCOtracking controller for a more challenging FCC unit with multiple inputs (Sect. 4.3).

Critical regions in mp-DO problem with first-order state path constraints
We consider the following problem: min u(t),x 1 (t),x 2 (t), 0≤t≤1 1 2 T 0.1242 0.0846 0.0846 0.9854 where the path constraint g(x) := 1.5x 1 + x 2 − 2 ≤ 0 is of first order with respect to dynamics. Moreover, the uncertainty is in the initial state conditions only, and the corresponding uncertainty domain is chosen as Two critical regions found on application of the algorithm in Sect. 3.2 are shown in Fig. 6(a), where part of the uncertainty domain is discarded due to infeasibility. The structure of the optimal control solutions in the critical region CR 1 consist of a unique interior arc, whereas those solutions in CR 2 are comprised of three arcs, an interior arc, followed by a boundary arc where the path constraint is binding, and a final interior arc constrained arc sharing the same control law as in CR 1 . Here, the boundary between CR 1 and CR 2 consists of those optimal control solutions where the path constraint activates at a single critical point along the time horizon [0, 1], and turns out to be nonlinear. For illustration, the optimal control and response trajectories along with the path constraint profile are shown in Fig. 6

Critical regions in mp-DO problem with second order state constraints
We now consider the following problem: where the path constraints are of second order. Here, uncertainty θ ∈ Θ := [0, 1] × [−1, 1] is present both in the initial conditions and the path/terminal constraints. A total of five critical regions is obtained for this mp-DO problem on application of the method presented in Sect. 3.2, as shown in Fig. 7(a). It turns out that this partition could in fact be extended to the entire right-half plane {θ 1 ≥ 0}. The critical region CR 1 corresponds to unconstrained optimal solutions. In CR 2 and CR 3 , the solutions are comprised of two interior arcs separated by a touch-and-go point for the path constraints x 1 (t) ≤ θ 1 and x 1 (t) ≥ −θ 1 , respectively, a behavior characteristic of higher-order path constraints; this behavior is illustrated for the optimal control solution for θ = (0.15, 0.8) ∈ CR 2 on Fig. 7(b), where the terminal constraint x 1 (1) ≤ 0 is also seen to be active. Finally, in CR 4 and CR 5 , the optimal solutions are comprised of three arcs, the first and last ones unconstrained and the middle one binding for the path constraints x 1 (t) ≤ θ 1 and x 1 (t) ≥ −θ 1 , respectively; this behavior is illustrated for the optimal control solution for θ = (0.1, 0.8) ∈ CR 4 on Fig. 7(c), where the path constraint x 1 (t) ≤ 0.1 remains active along the time interval [0.375, 0.625].

Multi-parametric NCO-tracking control of an FCC unit
This final case study considers a fluidized catalytic cracking (FCC) unit operated in partial combustion mode [41]. The objective is to steer the system to a given operating point, defined in terms of the mass fraction of coke on regenerated catalyst, C rc , and the regenerator dense bed temperature, T rg . The manipulated variables are the flow rate of air sent to the regenerator, F a , and the catalyst flow rate, F s . A linear input-output dynamic model is obtained via linearization and reduction of a first-principles nonlinear model around the equilibrium point C * rc = 5.207 × 10 −3 , T * rg = 965.4 K, T * ro = 776.9 K, T * cy = 988.1 K, and T * f = 400 K, where T f denotes the feed oil temperature. The control and state variables in this linear dynamic system are The optimization-based control problem is given by min u(t),x(t),y(t) with: The parameters θ in this mp-DO formulation correspond to uncertainty in C rc , T rg , and T f . We solve the mp-DO problem (29) with the approach presented in Sect  Fig. 8 both in a 3-d plot (Fig. 8(a)) and the corresponding 2-d projections for the parameter value θ 3 = 0 ( Fig. 8(b)); the discretized mp-QP counterparts for the latter will be presented later on (Fig. 10).
• The solutions in both CR 6 and CR 11 are comprised of three arcs, with a boundary arc where x 1 reaches its upper bound or its lower bound, respectively, located in between an initial interior arc and a final one. In both CR 5 and CR 10 , the solutions have the same constrained arc as CR 6 and CR 11 , respectively, yet without the final interior arc as the state constraint remains active until the terminal time.
• The solutions in CR 3 and CR 8 combine the previous two cases in CR 2 +CR 6 or CR 7 +CR 11 , and give rise to four arcs, starting with a boundary arc for u 2 (as in CR 2 or CR 7 ), followed by an interior arc, a boundary arc for x 1 (as in CR 6 or CR 11 ), and another interior arc; this case is illustrated in Fig. 9(c) for θ (3) = (6 × 10 −4 , 20, 0) ∈ CR 3 . The solutions in CR 4 and CR 9 present the same structure as CR 3 and CR 11 , respectively, but lack the final interior arc after the boundary arc where x 1 is active. For comparison, the critical regions shown in Fig. 10 are for the discretized mp-QP counterparts of (29) with either N = 20 or N = 50 time subintervals and for θ 3 = 0. In both cases, the central regions correspond to unconstrained solutions, like CR 1 in Fig. 8; all the other critical regions contain one or more active constraints along certain time subintervals, thereby approximating the six critical regions CR 2 -CR 11 in Fig. 8 and the nonlinear feedback control laws thereof in terms of affine control laws only. The actual switching structure gets better approximated as the time discretization is refined, but this also leads to a significant increase in the number of critical regions, namely 175 with N = 20 and 1687 with N = 50. In sum, the comparison with a classical mp-QP approach confirms the clear benefit offered by a continuoustime mp-DO approach in terms of a lesser number of critical regions, along with the ability to capture the underlying nonlinear feedback control nature. Closed-loop responses of the FCC unit based on the multi-parametric NCO-tracking controller derived from mp-DO problem (29) are shown in Fig. 11, starting from the uncertainty scenario (8 × 10 −4 , 15, 0) ∈ CR 6 at t = 0. A sampling period of ∆T = 10 −2 and signal-to-noise ratios of 50dB (Gaussian white noise) are considered in both cases, yet with different scenarios for the feed temperature. In the left plot ( Fig. 11(a)), a nominal feed temperature T /rm f is used, so θ 3 = 0 at all times; whereas the right plot ( Fig. 11(b) correspond to multiple step changes in the value of θ 3 as The resulting control and response trajectories (solid lines) present a good control performance compared to the nominal case (dashed lines: noise-free and infinite sampling).

Conclusions
This paper has introduced a framework for multi-parametric NCO tracking that exploits the multiparametric solution structure of an uncertain optimal control problem, without the need for applying a time discretization. In the special case of linear-quadratic optimal control problems, an algorithm has been (a) θ 3 = 0 all the time and 50dB signal-to-noise ratio (b) θ 3 with step changes and 50dB signal-to-noise ratio Figure 11: Closed-loop response of the multi-parametric NCO-tracking controller for problem (29).
proposed for characterizing the corresponding multi-parametric solution structure in terms of the exact critical regions and nonlinear feedback control laws. In practice, these feedback laws can be applied in a receding horizon manner, effectively resulting in a closed-loop, multi-parametric NCO-tracking controller for the system. In comparison to classical NCO tracking, this approach no longer requires the assumption of an invariant active set in the presence of uncertainty and extends the scope of NCO tracking to receding horizon control; whereas addressing the uncertain optimal control problem in its native, continuous-time form may lead to a dramatic reduction in the number of critical region compared to the classical mp-QP approach based on time discretization, due to the ability to capture the underlying nonlinear feedback control nature. These properties have been illustrated with several examples throughout the paper, including the two-input control problem of an FCC unit. Future work will focus on an extension of this approach to addressing (certain classes of) nonlinear systems. and an expression of the feedback control u(t) is given by if t ≤ t 1 (θ) −2 exp(−18.396[t − t 1 (θ)]), otherwise An implicit characterization of this critical region is obtained from (26) as CR (2) := θ ∈ [−2, 2] 2 −1 −1 θ ≤ −0.2440 .
• The point θ = (−1, −1) happens to be in the uncertainty domain Θ, yet outside CR (1) ∪ CR (2) . The structure S (3) of the optimal control at this point is also comprised of two arcs (N (3) t = 2), namely a boundary arc corresponding to an active upper input bound, followed by an interior arc. The determination of an explicit feedback control law and the characterization of CR (3) are similar to the previous case: • At this point, we have CR (1) ∪ CR (2) ∪ CR (3) = Θ and the procedure terminates with a total of three critical regions.