Multiparametric continuous-time optimal control via Pontryagin's Maximum Principle: explicit solutions and comparisons with discrete-time formulations

Multiparametric con tin uous-time optimal con trol via P on try agin’s Maxim um Principle: explicit solutions and comparisons with discrete-time form ulations ⋆ Lida Lamak ani a , b , Efstratios N. Pistik op oulos a , b a T exas A&M Ener gy Institute, T exas A&M University, Col le ge Station, USA b Artie McF errin Dep artment of Chemic al Engine ering, T exas A&M University, Col le ge Station, USA Abstract Mo del predictiv e con trol offers a pow erful framework for managing constrained systems, but its rep eated online optimization can b ecome computationally prohibitiv e. Multiparametric programming addresses this challenge b y precomputing optimal solutions offline, enabling real-time control through simple function ev aluation. While extensively developed for discrete-time systems, this approac h suffers from combinatorial gro wth in solution complexity as discretization is refined. This paper presen ts a systematic con tinuous-time m ultiparametric framew ork for linear-quadratic optimal control that directly solv es Pon tryagin’s optimalit y conditions without discretization artifacts. Through tw o illustrative examples, w e demonstrate that contin uous-time form ulations yield solutions with substan tially fewer critical regions than their discrete-time counterparts. Beyond this reduction in partition complexity , the contin uous-time approac h provides deeper insight in to system dynamics by explicitly identifying switc hing times and eliminating discretization artifacts that obscure the true structure of optimal con trol p olicies. Knowledge of the switching structure also accelerates online optimization methods b y providing analytical information ab out the solution top ology . Clear step-b y-step algorithms are provided for identifying switc hing structures, computing parametric switching times, and constructing critical regions, making the con tinuous-time framew ork accessible for practical implemen tation. Key wor ds: Multiparametric programming; Con tinuous-time optimal control; Explicit MPC; P ontry agin’s maximum principle; Critical regions; Switc hing times 1 In tro duction Mo del predictive control (MPC) is widely applied to complex constrained systems. The approach works b y using a mo del to predict how the system will b eha ve, then choosing control actions that optimize p erformance while sta ying within op erational limits [27]. This has made MPC esp ecially v aluable in settings where safety and physical constraints matter—chemical plants, auto- motiv e systems, rob otics, and finance [20,28,13,11,12]. MPC is t ypically implemen ted through a receding horizon strategy [27]. At each step, only the first con- trol action is applied before re-solving with the updated state, requiring a full optimization at ev ery sampling instan t. This computational burden naturally leads to the fol- ⋆ This pap er w as not presented at an y IF A C meet- ing. Corresp onding author E. N. Pistikopoulos E-mail. stratos@tam u.edu. Email addr esses: lida.lamakani@tamu.edu (Lida Lamak ani), stratos@tamu.edu (Efstratios N. Pistik op oulos). lo wing question: can the optimization b e solved once offline, rather than rep eatedly at each time step? Mul- tiparametric programming offers a solution by treat- ing the current state as a parameter, the optimization problem can b e solved for all feasible states offline. Pis- tik op oulos, Bemp orad, and colleagues [24,4] show ed that the optimal control law can b e expressed as a piecewise affine function ov er a partition of the state space. Online con trol then reduces to region identification and func- tion ev aluation, eliminating real-time optimization. This explicit MPC framework has b een extended to hybrid systems [6,21], uncertain and robust systems [3,30,18], and nonlinear systems [22,8,9], demonstrating substan- tial computational b enefits. The ma jority of explicit MPC researc h has consid- ered discrete-time formulations. These lead to finite- dimensional optimization problems that are computa- tionally tractable, explaining their p opularit y . A k ey lim- itation, how ever, is that most ph ysical systems are in- heren tly contin uous-time. Chemical reactors, vehicle dy- namics, rob otic manipulators, p ow er systems, and bi- ological processes are all naturally mo deled b y differ- en tial equations. Discretization in tro duces approxima- tions, and while finer time grids improv e accuracy , they also increase the problem dimension. Critically , the num- b er of critical regions in the parametric partition tends to grow substantially with the num b er of discretization p oin ts, increasing b oth offline computation and online lo okup cost. Can multiparametric metho ds b e developed directly for contin uous-time systems to av oid this com- binatorial explosion? A con tinuous-time formulation offers an alternative approac h that preserves the natural structure of the underlying dynamics. Rather than discretizing, the optimization is formulated directly in terms of the dif- feren tial equations. P on try agin’s Maxim um Principle [26] provides necessary conditions for optimality in the form of a t wo-point b oundary v alue problem. F or linear- quadratic systems, these conditions admit analytical solutions. In the presence of constraints, the optimal solution exhibits a piecewise structure: different analyt- ical expressions gov ern the solution ov er distinct time in terv als, with transitions occurring when constraints b ecome active or inactive. The switching times dep end on the initial state, which is precisely the parametric dep endence required for explicit MPC. When path constrain ts are active, the optimalit y con- ditions form a differential-algebraic equation (D AE) sys- tem combining the dynamics with algebraic constraint equations [7]. Biegler and colleagues ha ve developed sim ultaneous approaches that discretize DAE systems and solv e large-scale nonlinear programs [5]; sensitivity- based strategies hav e b een prop osed to reduce the asso- ciated online computational cost [34,35,10], though the discretization-dep enden t complexity discussed ab o ve remains. Bonvin and co work ers hav e developed comple- men tary metho ds for batch pro cesses that systemati- cally exploit the structure of optimalit y conditions to handle path constrain ts [32,31,14]. Their work demon- strates that understanding which constrain ts are active and when is crucial for both computational efficiency and physical insight—a p ersp ectiv e that aligns directly with the multiparametric contin uous-time framework dev elop ed here. Sakizlis and co w ork ers [30,29] pioneered this ap- proac h for linear con tin uous-time systems, demonstrat- ing that the parameter space can be partitioned in to critical regions where the switching times v ary contin u- ously as functions of the parameters while the sequence of active constraints remains fixed. Within each re- gion, the optimal tra jectories hav e explicit analytical forms—t ypically exponentials in time multiplied by functions of parameters—in contrast to the piecewise affine laws obtained in discrete-time formulations. Building on this foundation, the present work makes three contributions. First, we pro vide complete step- b y-step algorithms that transform the contin uous-time approac h from a theoretical concept into a practical computational metho d. W e reformulate the jump con- ditions using state contin uity instead of Hamiltonian con tin uit y , which is simpler to implement and more nu- merically stable. Second, we conduct the first system- atic study across multiple discretization levels, showing that contin uous-time form ulations consistently require far few er critical regions—often by an order of magni- tude—while discrete-time formulations exhibit com bi- natorial growth as discretization is refined. Third, w e demonstrate that contin uous-time solutions provide ex- plicit analytical expressions for switching times and op- timal tra jectories, offering direct ph ysical insigh t and enabling more efficient online optimization. In b oth examples, the contin uous-time partition holds at 5 critical regions regardless of accuracy requiremen ts, while the discrete-time coun t gro ws from 11 to o ver 60 as discretization is refined. This difference has direct practi- cal consequences: online region lo okup through 5 regions is substantially faster than through 60, which matte rs for systems op erating under tight computational budgets. The pap er is organized as follows. Section 2 formu- lates the contin uous-time optimal control problem and reviews optimality conditions. Section 3 develops the m ultiparametric framework. Section 4 presen ts tw o il- lustrativ e examples comparing contin uous-time and discrete-time approaches. Section 5 discusses implemen- tation and computational asp ects. Section 6 concludes. 2 Con tinuous-Time Optimal Control Structure This section form ulates both con tinuous-time and discrete-time optimal control problems and reviews the corresp onding optimality conditions. 2.1 Continuous-Time Optimal Contr ol Pr oblem F ollo wing the standard formulation in [29,25], the sys- tem dynamics are describ ed by ˙ x ( t ) = f  x ( t ) , u ( t )  , (1) where x ( t ) ∈ R n denotes the state vector and u ( t ) ∈ R n u is the con trol input. F or a given initial condition x ( t 0 ) = x , the con tin uous-time optimal con trol problem is p osed as min x ( · ) , u ( · ) J c ( x ) = Z t f t 0 ℓ  x ( t ) , u ( t )  dt + ϕ  x ( t f )  (2) s.t. ˙ x ( t ) = f  x ( t ) , u ( t )  , (3) g  x ( t ) , u ( t )  ≤ 0 , ∀ t ∈ [ t 0 , t f ] , (4) ψ  x ( t f )  ≤ 0 , (5) x ( t 0 ) = x. (6) F or comparison, w e also consider a discrete-time formu- lation obtained by sampling the dynamics with p erio d ∆ t , following standard discrete-time optimal co n trol and MPC formulations [22]. The sampled mo del is given by x k +1 = f d ( x k , u k ) , k = 0 , . . . , N − 1 , (7) where N = ( t f − t 0 ) / ∆ t and x 0 = x . A control horizon M ≤ N is in tro duced such that only the inputs { u k } M − 1 k =0 are treated as indep endent decision v ariables. F or k = M , . . . , N − 1, the control input is held constant, u k = u M − 1 , k = M , . . . , N − 1 . (8) 2 The corresp onding discrete-time optimal control prob- lem is min { x k ,u k } J d ( x ) = N − 1 X k =0 ℓ ( x k , u k ) + ϕ ( x N ) (9) s.t. x k +1 = f d ( x k , u k ) , k = 0 , . . . , M − 1 , (10) x k +1 = f d ( x k , u M − 1 ) , k = M , . . . , N − 1 , (11) g ( x k , u k ) ≤ 0 , k = 0 , . . . , M − 1 , (12) g ( x k , u M − 1 ) ≤ 0 , k = M , . . . , N − 1 , (13) ψ ( x N ) ≤ 0 , (14) 2.2 Assumptions Throughout the pap er, the following assumptions are adopted: • The system dynamics are linear and time-inv arian t. • State and input constraints are linear inequalities. • The p erformance index is quadratic, with Q ⪰ 0 and R ≻ 0. • A finite time horizon is considered. • The con trol horizon satisfies M ≤ N ; throughout this w ork M = N . 2.3 Pontryagin ’s Minimum Principle Conditions The necessary conditions for optimalit y of the con tin uous-time problem are given b y Pon tryagin’s Minim um Principle (PMP) [17,2]. The con tinuous-time optimal control problem is stated as min u ( t ) Φ( x ( t f )) + Z t f t 0 L ( x ( t ) , u ( t ) , t ) dt (15) sub ject to ˙ x ( t ) = f ( x ( t ) , u ( t ) , t ) , (16) g ℓ ( x ( t ) , u ( t ) , t ) ≤ 0 , ℓ = 1 , . . . , m, (17) ψ j ( x ( t f )) = 0 , j = 1 , . . . , p, x ( t 0 ) = x 0 . (18) The Hamiltonian function is defined as H ( x, u, λ, µ, t ) = L ( x, u, t ) + λ ⊤ ( t ) f ( x ( t ) , u ( t ) , t ) + m X ℓ =1 µ ℓ ( t ) g ℓ ( x, u, t ) . (19) The PMP conditions consist of the state equations ˙ x ( t ) = ∇ λ H = f ( x ( t ) , u ( t ) , t ) , (20) the costate equations ˙ λ ( t ) = −∇ x H , (21) and the optimality condition ∇ u H = 0 . (22) Constrain ts are enforced through complementary slac k- ness, µ ℓ ( t ) g ℓ ( x ( t ) , u ( t ) , t ) = 0 , ℓ = 1 , . . . , m, (23) together with feasibility conditions µ ℓ ( t ) ≥ 0 , g ℓ ( x ( t ) , u ( t ) , t ) ≤ 0 . (24) The terminal b oundary conditions for the costate v ari- ables are given by the transversalit y conditions asso ci- ated with terminal equality constraints [33], λ ( t f ) = ∇ x Φ  x ( t f )  + p X j =1 ν j ∇ x ψ j  x ( t f )  . (25) When the activit y status of constrain ts c hanges, the op- timal solution ma y exhibit corner p oints at which the classical W eierstrass–Erdmann conditions apply [15]. In particular, for problems with con tin uous dynamics and no switc hing costs, the costate v ariables and the Hamil- tonian are con tin uous across the switc hing time t s , lead- ing to the jump conditions λ ( t + s ) = λ ( t − s ) , H ( t + s ) = H ( t − s ) , (26) where t s ∈ ( t 0 , t f ) denotes a switching p oin t at which the activit y status of one or more path constraints changes. The structure and in terpretation of switching p oin ts are discussed in detail in the subsequent section. F or comparison, a discrete-time coun terpart of the ab o ve problem is obtained by time discretization ov er a finite horizon. Let t k = t 0 + k ∆ t , k = 0 , . . . , N , and define x k ≈ x ( t k ) and u k ≈ u ( t k ). The discrete-time op- timal control problem is written as min { x k ,u k } Φ( x N ) + N − 1 X k =0 L ( x k , u k , k ) (27) sub ject to x k +1 = f d ( x k , u k ) , k = 0 , . . . , N − 1 , (28) g ℓ,k ( x k , u k ) ≤ 0 , ℓ = 1 , . . . , m, k = 0 , . . . , N − 1 , (29) ψ j ( x N ) = 0 , j = 1 , . . . , p, x 0 = x 0 . (30) Unlik e the con tinuous-time case, the discrete-time prob- lem is finite-dimensional and its optimalit y conditions can b e written using the Karush–Kuhn–T uc ker (KKT) 3 conditions. The asso ciated Lagrangian function is L = Φ( x N ) + N − 1 X k =0 L ( x k , u k , k ) + N − 1 X k =0 λ ⊤ k +1  x k +1 − f d ( x k , u k )  + N − 1 X k =0 m X ℓ =1 µ ℓ,k g ℓ,k ( x k , u k ) + p X j =1 ν j ψ j ( x N ) . (31) The KKT conditions consist of the stationarity condi- tions ∇ x k L = 0 , k = 1 , . . . , N , (32) ∇ u k L = 0 , k = 0 , . . . , M − 1 , (33) ∂ L ∂ ν j = 0 , j = 1 , . . . , p, (34) together with feasibility of the constraints x k +1 − f d ( x k , u k ) = 0 , k = 0 , . . . , N − 1 , (35) g ℓ,k ( x k , u k ) ≤ 0 , ℓ = 1 , . . . , m, k = 0 , . . . , N − 1 , (36) ψ j ( x N ) = 0 , j = 1 , . . . , p, (37) dual feasibility µ ℓ,k ≥ 0 , ℓ = 1 , . . . , m, k = 0 , . . . , N − 1 , (38) and complementary slackness µ ℓ,k g ℓ,k ( x k , u k ) = 0 , ℓ = 1 , . . . , m, k = 0 , . . . , N − 1 . (39) A k ey distinction betw een the tw o form ulations lies in the nature of their optimalit y conditions. In the con tin uous-time case, Pon tryagin’s Minimum Principle yields a system of ordinary differential equations with costates and multipliers ev olving contin uously ov er time. In contrast, the discrete-time formulation leads to a fi- nite set of algebraic optimality conditions through the KKT framew ork. This fundamental difference b etw een differen tial and algebraic conditions directly influences the structure of switching times and critical regions in con tin uous- and discrete-time formulations. 2.4 Switching Points and Jump Conditions In contin uous-time constrained optimal control, the optimal solution typically changes its structure as con- strain ts become activ e or inactiv e o ver time. W e refer to time in terv als where the same constrain ts remain active as ar cs . During each arc, the state and control tra jecto- ries follo w a sp ecific analytical form dictated by which constrain ts are curren tly enforced. When a constraint switc hes from active to inactive (or vice versa), the so- lution structure c hanges; these transition p oints in time are called switching p oints or c orners . F or constrained arcs, it is useful to distinguish b e- t w een the time at which a constraint b ecomes activ e and the time at which it is released. Each constrained arc is therefore characterized by tw o switching instants: Fig. 1. Schematic illustration of switching points and con- strained arcs. The times t st and t sx denote the entry and exit points of a constrained arc, resp ectively . Changes in the active constrain t set are accompanied b y corresp onding c hanges in the asso ciated Lagrange m ultipliers. T able 1 Constrain t signs and Lagrange m ultipliers at en try and exit switc hing p oints. Corner Constraint signs Lagrange multipliers En try/exit of constraint g 1 t − 1 t g 1 < 0 , g 2 < 0 µ 1 = 0 , µ 2 = 0 t + 1 t g 1 = 0 , g 2 < 0 µ 1 > 0 , µ 2 = 0 t − 1 x g 1 = 0 , g 2 < 0 µ 1 > 0 , µ 2 = 0 t + 1 x g 1 < 0 , g 2 < 0 µ 1 = 0 , µ 2 = 0 En try/exit of constraint g 2 t − 2 t g 1 < 0 , g 2 < 0 µ 1 = 0 , µ 2 = 0 t + 2 t g 2 = 0 , g 1 < 0 µ 1 = 0 , µ 2 > 0 t − 2 x g 2 = 0 , g 1 < 0 µ 1 = 0 , µ 2 > 0 t + 2 x g 1 < 0 , g 2 < 0 µ 1 = 0 , µ 2 = 0 • Entry p oint ( t st ): the time at which the s -th con- strain t b ecomes active. • Exit p oin t ( t sx ): the time at which the s -th con- strain t is deactiv ated. Figure 1 illustrates these concepts schematically , while T able 1 summarizes the relationship b etw een switc hing p oin ts, constraint activit y , and the asso ciated Lagrange m ultipliers. F or simplicit y , when there is only one switching point w e denote it by t s rather than t st or t sx . F rom the standp oin t of the P ontry agin Minim um Principle, switching p oints require sp ecial treatmen t to ensure consistency of the optimality conditions across adjacen t arcs. F or problems with contin uous dynamics and order-zero path constrain ts, t wo mathematically equiv alen t formulations can b e emplo y ed. The classical approac h enforces contin uity of both the costate and the Hamiltonian, λ ( t + s ) = λ ( t − s ) , H ( t + s ) = H ( t − s ) , (40) while an alternative approac h enforces contin uity of the state and the costate, x ( t + s ) = x ( t − s ) , λ ( t + s ) = λ ( t − s ) . (41) Both formulations provide tw o scalar conditions per switc hing p oint and lead to identical optimal solutions. 4 Ho w ever, the state-based form ulation (41) offers practi- cal adv antages in numerical implementations, as state v ariables are directly integrated from the dynamics and exhibit smo other b ehavior near switching p oin ts com- pared to the Hamiltonian, which inv olves comp osite expressions of states, costates, and multipliers. The following prop osition formalizes the connection b et ween these tw o formulations and establishes condi- tions under which the state-based approach is sufficient. Prop osition 2.1 (State-matc hing c haracterization of switc hing p oints) . Consider a c ontinuous-time optimal c ontr ol pr oblem with p ath c onstr aints of or der zer o (de- p ending only on x , u , and t ). Assume c ontinuous dy- namics without impulsive terms and no switching c osts. Then, for a pr escrib e d ar c se quenc e, the switching times c an b e determine d by enfor cing the jump c onditions x ( t − s ) = x ( t + s ) , λ ( t − s ) = λ ( t + s ) , (42) to gether with the c onstr aint activity c onditions. Under these assumptions, Hamiltonian c ontinuity fol lows auto- matic al ly and do es not ne e d to b e imp ose d as an addi- tional e quation. Pr o of. State contin uity follows directly from the ab- sence of impulsive dynamics. Costate contin uity follows from the W eierstrass–Erdmann corner conditions, which guaran tee contin uous adjoint v ariables across switch- ing p oints for order-zero constraints without switching costs. A t the switching instant, the constraint b oundary is reac hed: g i ( x s , u, t s ) = 0. Since x s and λ s are identi- cal from both sides, the stationarit y condition ∂ H /∂ u = 0 together with the constraint equation determine a unique con trol u s under standard regularit y assumptions (strict conv exity of H in u and linear independence of activ e constraint gradients). Thus, u ( t − s ) = u ( t + s ) = u s . The Hamiltonian is H = L + λ ⊤ f + µ ⊤ g . At the switch- ing p oint, g ( x s , u s , t s ) = 0, so the µ ⊤ g term v anishes regardless of µ . Since x s , u s , and λ s are all contin uous across t s , we hav e H ( t − s ) = H ( t + s ). R emark 2.1 . An equiv alent alternative formulation uses Hamiltonian contin uity instead of state contin uity: H ( t − s ) = H ( t + s ) , λ ( t − s ) = λ ( t + s ) . (43) Both formulations provide tw o scalar conditions per switc hing p oint. This work adopts the state-based for- m ulation (42) b ecause the state v ariables x ha v e direct ph ysical meaning and yield more n umerically stable computations, particularly for autonomous systems where the Hamiltonian is constan t along the optimal tra jectory . When the assumptions of Prop osition 2.1 are violated, mo dified jump conditions apply: • Higher-order constraints (e.g., g ( ˙ x ) or g ( ¨ x )): The costate exhibits a jump discon tin uit y giv en b y ∆ λ ( t s ) = − µ ( t s ) ∂ g /∂ ˙ x . • Impulsive dynamics : The state exhibits a jump discon tin uity , and the Hamiltonian-based formulation m ust b e used instead of state matching. • Switching costs : The Hamiltonian exhibits a jump discon tin uity given by ∆ H ( t s ) = − ∂ c/∂ t s , where c denotes the switching cost function. 3 Multiparametric Con tinuous-Time Optimal Con trol The contin uous-time optimal con trol problem is for- m ulated in a multiparametric framework where param- eters represent either initial states or data uncertaint y within a given set. W e fo cus on initial state parametriza- tion, c haracterizing ho w optimal con trol la ws, switching times, and arc structures v ary with parameters. Rather than solving individual parameter v alues re- p eatedly , the multiparametric approach treats the en- tire problem family in a unified framework. Each ad- missible parameter v alue yields an optimal solution sat- isfying Pon tryagin’s optimalit y conditions. W e aim to understand ho w solutions ev olv e with parameter v aria- tions and iden tify parameter space regions where solu- tion structure remains qualitatively unchanged. In contin uous-time systems, parameter v ariations di- rectly affect the active constraint sequence and switching times, creating piecewise optimal solutions whose struc- ture dep ends on parameter v alues. By analyzing this de- p endence explicitly , w e can express optimal control laws and state tra jectories as explicit parameter functions o v er paramete r subspace regions. This dev elopment extends classical contin uous-time optimal con trol theory through the m ultiparametric framew ork, pro viding a foundation for critical region analysis and explicit solution representations. 3.1 Par ametric F ormulation In the multiparametric setting, we treat the contin uous- time optimal control problem as a family indexed by parameter v ector θ , studying ho w optimal solutions v ary across an admissible set rather than solving for fixed parameter v alues. Throughout this work, θ represen ts the initial state of the system, denoted as x ( t 0 ) = θ ≡ x 0 , θ ∈ Θ ⊂ R n , (44) where Θ denotes a compact set of admissible initial con- ditions. The parametric problem then takes the form min x ( · ) , u ( · ) J ( x, u ; θ ) = Z t f t 0 ℓ  x ( t ) , u ( t )  dt + ϕ  x ( t f )  (45) s.t. ˙ x ( t ) = f  x ( t ) , u ( t )  , (46) g  x ( t ) , u ( t )  ≤ 0 , ∀ t ∈ [ t 0 , t f ] , (47) ψ  x ( t f )  ≤ 0 , (48) x ( t 0 ) = θ , θ ∈ Θ . (49) F or each θ ∈ Θ, this formulation yields a standard con tin uous-time optimal con trol problem satisfying 5 P on try agin’s optimality conditions. The multiparamet- ric framework simply provides a unified wa y to represent the entire family of problems across differen t parameter v alues without changing the underlying dynamics or constrain ts. 3.2 Par ametric Switching Structur es W e no w examine how the structure of the optimal solution, particularly the switching instants, dep ends on the initial state θ = x 0 . All inequality constrain ts are written in the compact form g ( x ( t ) , u ( t )) ≤ 0 , (50) with asso ciated Lagrange multipliers satisfying µ ( t ) ≥ 0 . (51) Switc hing o ccurs exclusively due to changes in con- strain t activity . A switching instant is defined as a time at which the activ e set of constrain ts changes, or equiv a- len tly when at least one Lagrange multiplier transitions b et ween zero and a p ositive v alue. The time horizon can therefore b e decomp osed as t 0 < t 1 < · · · < t N s < t f , (52) where the num b er of switching instants N s dep ends on the initial state x 0 . Each switching instant is treated as a function of the initial state, t s = t s ( x 0 ) , s = 1 , . . . , N s . (53) A t a switc hing instant t s , the jump conditions for path-constrained problems without impulsive effects re- quire contin uity of the state and adjoint v ariables, x ( t + s ) = x ( t − s ) , (54) λ ( t + s ) = λ ( t − s ) . (55) The junction conditions at a switching instant are ex- pressed through constrain t activit y and the correspond- ing multiplier b ehavior. En try of the constrain t asso ci- ated with switching even t s is characterized by g ( t − s ) < 0 , µ ( t − s ) = 0 , g ( t + s ) = 0 , µ ( t + s ) > 0 , (56) whereas exit of the constraint is characterized by g ( t − s ) = 0 , µ ( t − s ) > 0 , g ( t + s ) < 0 , µ ( t + s ) = 0 . (57) All remaining constrain ts satisfy feasibilit y and comple- men tary slac kness on both sides of the switching instant. A parametric switc hing structure is defined as a fixed ordered sequence of active constraint sets ov er the time horizon. A critical region R ⊂ Θ is a subset of the initial- state space such that for ev ery x 0 ∈ R the same switch- ing structure is observed and the jump conditions (54)– (55) and junction conditions (56)–(57) remain v alid. Algorithm 1 Identification and fitting of parametric switc hing times Input: Switching structure S identified at a seed initial state x ⋆ 0 , and sampled initial states { x ( r ) 0 } N S r =1 restricted to the domain where S is observed Output: Empirical switc hing-time functions { ˆ t s ( x 0 ) } N s s =1 asso ciated with S Indexing: r indexes sampled initial states; s indexes switc hing even ts (1) Solve the con tinuous-time optimality conditions at x ⋆ 0 to iden tify the switching structure S and the ordered switching even ts s = 1 , . . . , N s . (2) F or each sample r = 1 , . . . , N S : (a) Solve the optimality system consisten t with S at x ( r ) 0 . (b) F or eac h switching even t s = 1 , . . . , N s , com- pute the exact switching time t ( r ) s = t s ( x ( r ) 0 ) b y enforcing the jump conditions (54)–(55) and the junction conditions (56) or (57). (c) Retain t ( r ) s if all referenced conditions are sat- isfied. (3) F or each switching ev en t s = 1 , . . . , N s , fit an em- pirical function ˆ t s ( x 0 ) using the retained sample pairs { ( x ( r ) 0 , t ( r ) s ) } . (4) return { ˆ t s ( x 0 ) } N s s =1 . Since the num b er and lo cations of switching instants are not kno wn beforehand, a structure-driven pro cedure is adopted. F or a fixed switching structure, exact switch- ing times t s ( x 0 ) are computed for selected initial states b y solving the optimalit y system together with the jump and junction conditions. These exact v alues are then used to construct empirical approximations ˆ t s ( x 0 ) de- scribing the dep endence of switching instants on the ini- tial state within the same structure. The pro cedure is summarized in Algorithm 1. Lemma 3.1 (Existence, uniqueness, and contin uity of switc hing time) . F or a fixe d switching structur e S with a single switching event, ther e exists a unique switching time t s ( x 0 ) ∈ ( t 0 , t f ) for e ach x 0 ∈ R , and the switching time function t s : R → [ t 0 , t f ] is c ontinuous over the interior of the critic al r e gion R . Pr o of. Step 1: Existence. F or each x 0 ∈ R , the switch- ing time t s ( x 0 ) is defined as the time when the Lagrange m ultiplier reaches zero: µ ( t s , x 0 ) = 0 . (58) F or linear-quadratic problems, the multiplier µ ( t, x 0 ) is obtained by solving the coupled system of linear ordinary differen tial equations gov erning the state and costate, with initial conditions determined by x 0 . Since the crit- ical region R is defined by a fixed switc hing structure, a switc hing time exists for each x 0 ∈ R b y construction. Step 2: F orm of the multiplier. F or linear systems, the state and costate ev olv e according to the Hamilto- nian system, whic h is a 2 n -dimensional linear ordinary 6 differen tial equation. The m ultiplier µ ( t, x 0 ) is deter- mined algebraically from the stationarity condition and can therefore b e expressed as a linear combination of the fundamen tal solutions: µ ( t, x 0 ) = 2 n X i =1 α i ( x 0 ) e σ i t , (59) where σ i are the 2 n eigenv alues of the Hamiltonian sys- tem matrix and α i ( x 0 ) are coefficients that dep end con- tin uously on the initial state x 0 . This exp onential form follo ws directly from the linearit y of the state-costate system. Step 3: Uniqueness within the critical region. Within a critical region R , the switching structure is fixed by definition. F or a single switching ev en t, the con- strain t transitions from activ e to inactiv e (or vice v ersa) exactly once. Supp ose µ ( t, x 0 ) = 0 had m ultiple solutions in ( t 0 , t f ). Then the m ultiplier w ould cross zero multiple times, im- plying m ultiple constrain t activ ations or deactiv ations. This w ould con tradict the assumption of a fixed switch- ing structure with a single switching even t. Therefore, for each x 0 ∈ R , there exists a unique t s ( x 0 ) ∈ ( t 0 , t f ) satisfying µ ( t s ( x 0 ) , x 0 ) = 0. Step 4: Non-v anishing time deriv ative. At the switc hing time t s , the constrain t transitions b etw een ac- tiv e and inactive states. W e show that ∂ µ ∂ t ( t s , x 0 )  = 0. F rom Step 2, we hav e µ ( t, x 0 ) = P 2 n i =1 α i ( x 0 ) e σ i t , whic h is smooth and has a w ell-defined deriv ative at ev- ery p oint. Supp ose, for the sake of con tradiction, that ∂ µ ∂ t ( t s , x 0 ) = 0. Since µ ( t s , x 0 ) = 0 by definition of the switc hing time, this gives: µ ( t s ) = 0 and dµ dt ( t s ) = 0 . (60) Consider the T aylor expansion of µ around t s : µ ( t ) = µ ( t s )+ dµ dt ( t s )( t − t s )+ O  ( t − t s ) 2  = O  ( t − t s ) 2  . (61) This implies that µ ( t ) remains near zero with second- order smallness in a neighborho o d of t s . Ho w ev er, within the critical region R , the switching structure requires that on one side of t s the constraint is active (with µ > 0), and on the other side it is inactiv e (with µ = 0 by complemen tarit y). If dµ dt ( t s ) = 0, then µ cannot transition immedi- ately from a p ositiv e v alue to zero. It w ould instead ho v er near zero, requiring higher-order deriv atives to v anish as w ell. But from the exp onen tial form µ ( t, x 0 ) = P 2 n i =1 α i ( x 0 ) e σ i t , if µ and all its deriv atives v anish at t s , then µ ≡ 0 for all t , contradicting the existence of an active constraint phase. Therefore, ∂ µ ∂ t ( t s , x 0 )  = 0. Step 5: Con tinuit y . Since µ ( t, x 0 ) is con tin uously differen tiable in both argumen ts and ∂ µ ∂ t ( t s , x 0 )  = 0, the Implicit F unction Theorem [19] guarantees that there exists a unique con tin uous function t s ( x 0 ) in a neigh- b orho o d of any x ∗ 0 ∈ in t( R ) satisfying µ ( t s ( x 0 ) , x 0 ) = 0. Since this holds for every p oint in the interior of R , the switc hing time function t s : R → [ t 0 , t f ] is con tin uous o v er int( R ). R emark 3.1 . The contin uity established in Lemma 3.1 has imp ortant practical implications for the represen- tation of switching-time functions. By the W eierstrass Appro ximation Theorem, an y contin uous function on a compact set can b e uniformly approximated b y p olyno- mials to arbitrary precision. Since critical regions are b ounded subsets of the parameter space and t s ( x 0 ) is con tin uous ov er eac h region, p olynomial approxima- tions are theoretically guaranteed to con v erge to the true switching-time function. When an analytical expression is av ailable, it should b e used directly . When analytical deriv ation is in- tractable, the W eierstrass theorem justifies the em- pirical fitting pro cedure in Algorithm 1. In practice, lo w-degree p olynomials achiev e high accuracy ov er the relativ ely narro w parameter ranges of individual critical regions. Alternative approximation schemes—suc h as rational functions, splines, or other basis functions— ma y b e emplo y ed dep ending on the complexity of the switc hing-time dep endence and the desired balance b et ween accuracy and computational simplicity . 3.3 Critic al R e gions and Explicit Solutions The parametric switching structures identified in the previous subsection induce a partition of the initial-state space in to subsets within which the qualitative structure of the optimal solution remains unchanged. These sub- sets are referred to as critic al r e gions and provide the basis for constructing explicit parametric solutions. F or a fixed switching structure, a critical region is the set of initial states x 0 ∈ Θ for whic h the same constraints remain active and inactive o ver the entire time horizon. In practice, this is enforced by t w o simple criteria that mirror the discrete-time definition of critical regions (for example, in multiparametric quadratic programming or explicit MPC), with the only difference b eing that the conditions must hold ov er a contin uous-time horizon: • Inactive constrain ts remain inactive. • Active Lagrange multipliers remain p osi tive. Accordingly , a critical region R is characterized by parametric inequalities ev aluated along the optimal tra- jectory . Let ˜ g ( · ) denote the vector of inactive inequalit y constrain ts, and let ˜ µ ( · ) and ˜ ν ( · ) denote the Lagrange m ultipliers asso ciated with active constraints and their en try points, resp ectively . The interior of a critical re- gion is defined by ˜ g  x ( t, t s ( x 0 ) , x 0 ) , u ( t, t s ( x 0 ) , x 0 ) , x 0  < 0 , (62) ˜ µ  t, t s ( x 0 ) , x 0  > 0 , (63) ˜ ν  t s ( x 0 ) , x 0  > 0 . (64) The inequalities (62)–(64) are time-dep endent. How- ev er, to obtain an explicit description of a critical re- gion in the initial-state space, it is sufficient to enforce 7 them at their most restrictive p oints in time. This yields explicit parametric expressions that define the region b oundaries. F or each inactive constraint i , the b oundary is given b y ¯ G i ( x 0 ) = max t ∈ [ t 0 ,t f ] ˜ g i  x, u, x 0 ) , x 0  , (65) and the condition ¯ G i ( x 0 ) < 0 defines the interior of the region with resp ect to that constraint. F or constraints that are active o v er at least one con- strained arc [ t − i,s ( x 0 ) , t + i,s ( x 0 )], the inactiv e p ortions of the horizon must also satisfy ¯ G i ( x 0 ) = max t ∈ [ t 0 ,t f ] t / ∈ [ t − i,s ( x 0 ) ,t + i,s ( x 0 )] ˜ g i  x, u, x 0  , (66) Finally , for each active constraint i , p ositivity of the asso ciated Lagrange multiplier at all en try p oin ts must b e preserved, which yields the b oundary condition ¯ µ i ( x 0 ) = min t ∈ [ t − i,s ( x 0 ) ,t + i,s ( x 0 )] µ i  x, u, x 0  , (67) Equations (65)–(67) provide explicit parametric ex- pressions for the b oundaries of the critical regions. Cross- ing any of these b oundaries implies that at least one pre- viously inactive constraint b ecomes active (or an active m ultiplier loses p ositivity), which c hanges the switc hing structure and therefore leads to a differen t explicit solu- tion. Within a critical region, the switc hing structure is fixed and the optimality system admits a uniform rep- resen tation. As a result, the optimal control, state, and adjoin t tra jectories can b e expressed explicitly as func- tions of time and the initial state. The optimal con trol la w takes the piecewise form u ∗ ( t, x 0 ) =          u 1 ( t, x 0 ) , t ∈ [ t 0 , t 1 ( x 0 )) , u 2 ( t, x 0 ) , t ∈ [ t 1 ( x 0 ) , t 2 ( x 0 )) , . . . u N s +1 ( t, x 0 ) , t ∈ [ t N s ( x 0 ) , t f ] , (68) where each arc-wise expression is obtained from the Hamiltonian minimization condition under the fixed activ e set asso ciated with the region. Figure 2 shows how the tw o algorithms in teract: Al- gorithm 2 drives the global exploration, calling Algo- rithm 1 lo cally to fit switching times for each identified structure. 4 Illustrativ e Examples This section presents t wo illustrative examples to demonstrate the formulation and solution of m ulti- parametric dynamic optimization problems in b oth con tin uous-time (CT) and discrete-time (DT) settings. F or each example, the contin uous-time problem is first solv ed analytically , leading to an explicit partition of Algorithm 2 Solution of the multiparametric dynamic optimization problem Input: Parameter domain Θ, initial parameter seed θ (0) ∈ Θ Output: Set of critical regions and exact parametric so- lutions (1) Initialize the unexplored parameter set Θ rest ← Θ. (2) While Θ rest  = ∅ do : (a) Select a fixed parameter p oint θ ( i ) ∈ Θ rest . (b) Solve the dynamic optimization problem at θ ( i ) to identify the optimal arc structure, active constrain ts, and switching p oints. (c) Derive exact analytical expressions for the state, adjoin t, control, and multipliers corre- sp onding to this structure. (d) Construct symbolic b oundaries of the critical regions b y enforcing feasibility , complementar- it y , and multiplier sign conditions using (65)– (67). (e) Define the asso ciated critical region R i ⊂ Θ as the set of parameters satisfying these condi- tions. (f ) Remov e the iden tified region from the unex- plored set: Θ rest ← Θ rest \ R i . (3) return the collection of critical regions and their exact parametric solutions. Parameter domain Θ Θ rest ← Θ Alg. 2 Global region exploration Select seed θ ( i ) ∈ Θ rest Solve DO at θ ( i ) Identify switching structure S Alg. 1 Fit switching times ˆ t s ( x 0 ) Derive explicit arc-wise expressions Construct critical region R i Θ rest ← Θ rest \ R i Return {R i } and explicit solutions if Θ rest = ∅ Fig. 2. Interaction b etw een Alg. 2 (global exploration) and Alg. 1 (lo cal switching-time fitting). the parameter space in to critical regions with region- sp ecific optimal tra jectories. Once the critical region corresp onding to a given initial condition is identified, 8 the optimal state and control tra jectories are obtained b y direct ev aluation of the asso ciated expressions, with switc hing handled through region-dep endent switching times when applicable. All contin uous-time computations are p erformed sym- b olically using SymPy , with switching-time expressions ev aluated via lambdify . When closed-form expressions are not retained, sampled solutions are p ost-pro cessed using NumPy to construct low-degree p olynomial approx- imations for rep orting. The discrete-time coun terparts are formulated as multiparametric quadratic programs and solv ed using the PPOPT to olb ox [16]. Sp ecifically , w e employ the combinatorial parallel algorithm imple- men ted in PPOPT , whic h systematically explores the space of activ e sets to obtain the corresp onding critical- region partitions and affine control laws. 4.1 Example 1: One-State Inte gr ating System with Mixe d Output and Input Constr aints This example illustrates the form ulation and solution of a multiparametric dynamic optimization problem for a simple contin uous-time system with path constraints. 4.1.1 Pr oblem Statement min u ( · ) 1 2 x ( T ) 2 + 1 2 Z T 0  x ( t ) 2 + u ( t ) 2  dt s.t. ˙ x ( t ) = − u ( t ) , t ∈ [0 , T ] , y ( t ) = x ( t ) + u ( t ) , t ∈ [0 , T ] , − 1 ≤ y ( t ) ≤ 1 , t ∈ [0 , T ] , u ( t ) ≤ 2 , t ∈ [0 , T ] , x (0) = x 0 , − 2 ≤ x 0 ≤ 2 , T = 2 . (69) The initial state x 0 is the parameter, and the aim is to deriv e an explicit representation of the optimal control la w ov er Θ. The construction is presented b elow. 4.1.2 Continuous-Time Multip ar ametric Solution The initial state x 0 is treated as the scalar parameter on the domain Θ := { x 0 ∈ R | − 2 ≤ x 0 ≤ 2 } . The ob jective is to deriv e an explicit parametric descrip- tion of the optimal tra jectories and their switching struc- ture ov er Θ. Iter ation 1. (1) Initialize the unexplored parameter set: Θ rest ← Θ . (2a) Fix a seed parameter in the unexplored set: x (0) 0 := 0 ∈ Θ rest . (2b) At the seed p oin t x (0) 0 = 0, the optimality condi- tions yield a tra jectory that satisfies all path and input constraints throughout [0 , 2], so no switc h- ing p oints arise and the solution is a s ingle uncon- strained arc. (2c) Derive the exact analytical expressions for the state, con trol, adjoin t, and constraint-related quantities corresp onding to this arc. The Hamiltonian is H = 1 2 ( x 2 + u 2 ) + λ ( − u ) , whic h gives the stationarity condition u ( t ) = λ ( t ). The state and adjoint equations are ˙ x ( t ) = − u ( t ) , ˙ λ ( t ) = − x ( t ) , with terminal condition λ (2) = x (2). Solving yields x ( t ) = x 0 e − t , u ( t ) = x 0 e − t , λ ( t ) = x 0 e − t , t ∈ [0 , 2] . The path constraint function is y ( t ) = x ( t ) + u ( t ) , and since the constraint is inactive on this arc, the corresp onding multiplier satisfies µ ( t ) = 0 ∀ t ∈ [0 , 2] . (2d) Enforce feasibilit y conditions sym b olically o ver the horizon. The output is y ( t ) = x ( t ) + u ( t ) = 2 x 0 e − t . Requiring − 1 ≤ y ( t ) ≤ 1 for all t ∈ [0 , 2] gives − 1 ≤ 2 x 0 e − t ≤ 1 . Since e − t ∈ (0 , 1], feasibility reduces to − 1 ≤ 2 x 0 ≤ 1 ⇐ ⇒ − 0 . 5 ≤ x 0 ≤ 0 . 5 . (2e) Define the asso ciated critical region: CR03 CT = { x 0 ∈ Θ | − 0 . 5 ≤ x 0 ≤ 0 . 5 } , with a single unconstrained arc on [0 , 2]. (2f ) Up date the unexplored parameter set: Θ rest ← [ − 2 , − 0 . 5) ∪ (0 . 5 , 2] . (3) Since Θ rest  = ∅ , pro ce ed to the next iteration. Iter ation 2. (2a) Fix a new seed parameter in the unexplored set: x (1) 0 := − 0 . 8 ∈ Θ rest . 9 (2b) At x (1) 0 = − 0 . 8, the optimality conditions show that the unconstrained solution violates the lo wer out- put b ound at t = 0. The optimal tra jectory there- fore has tw o arcs: an initial constrained arc with the lo w er b ound active, follow ed by an unconstrained arc after the constrain t is released at a single switch- ing p oint. (2c) Derive exact analytical expressions for the state, con trol, and constraint-related quan tities for the iden tified tw o-arc structure. Define the low er out- put constraint as g ( t ) = − x ( t ) − u ( t ) − 1 ≤ 0 . On the first arc the constraint is active, hence g ( t ) = 0 , µ ( t ) ≥ 0 . whic h implies u ( t ) = − 1 − x ( t ) , ˙ x ( t ) = 1 + x ( t ) . With x (0) = x 0 , the constrained-arc solution is x 0 ( t ) = ( x 0 +1) e t − 1 , u 0 ( t ) = − ( x 0 +1) e t , t ∈ [0 , t s ] . On the second arc the constraint is inactive, so that g ( t ) ≤ 0 , µ ( t ) = 0 . The unconstrained optimal feedback applies, u ( t ) = x ( t ) , ˙ x ( t ) = − x ( t ) , yielding, by contin uity at t s , x 1 ( t ) = x ( t s ) e − ( t − t s ) , u 1 ( t ) = x ( t s ) e − ( t − t s ) , t ∈ [ t s , 2] . (2d) Enforce the symbolic b oundary and junction condi- tions. The switc hing time t s corresp onds to lea ving the active constraint. A t this p oin t the tra jectory m ust satisfy the constraint b oundary and b e com- patible with the unconstrained feedback law. On the unconstrained arc, u ( t ) = x ( t ) implies g ( t ) = − x ( t ) − u ( t ) − 1 = − 2 x ( t ) − 1 . Enforcing g ( t s ) = 0 therefore yields − 2 x ( t s ) − 1 = 0 = ⇒ x ( t s ) = − 1 2 . Substituting the constrained-arc expression gives ( x 0 +1) e t s − 1 = − 1 2 = ⇒ t s ( x 0 ) = ln  1 2( x 0 + 1)  . Requiring t s ( x 0 ) ∈ [0 , 2] yields − 1 + 1 2 e 2 ≤ x 0 ≤ − 1 2 . T able 2 Activ e-arc structures of the contin uous-time critical regions and their corresp onding parameter interv als. Region x 0 in terv al Arc structure CR01 CT [ − 1 . 2707 , − 0 . 9323] y min activ e (full horizon) CR02 CT [ − 0 . 9323 , − 0 . 5] y min activ e → unconstrained CR03 CT [ − 0 . 5 , 0 . 5] unconstrained (full horizon) CR04 CT [0 . 5000 , 0 . 9323] y max activ e → unconstrained CR05 CT [0 . 9323 , 2] y max activ e (full horizon) (2e) Define the asso ciated critical region: CR02 CT =  x 0 ∈ Θ     − 1 + 1 2 e 2 ≤ x 0 ≤ − 0 . 5  . (2f ) Remov e this region from the unexplored set. The remaining p ortion of the negative side is [ − 2 , − 0 . 5] \ CR02 CT = [ − 2 , − 1 + 1 2 e 2 ) . (3) At this stage, tw o critical regions hav e b een con- structed explicitly: CR03 CT and CR02 CT . By symmetry , the positive half of the parameter space yields analogous regions with opp osite signs. The c om- plete pro cedure identifies five critical regions. The activ e-arc structures associated with the contin uous- time critical regions are summarized in T able 2. These structures iden tify whic h constraints are activ e o v er the time horizon for eac h parameter interv al. Giv en the active-arc structure of eac h region, the corresp onding optimal tra jectories can b e deriv ed in closed form. The resulting closed-form expressions for the optimal state and control tra jectories are rep orted in T able 3. Switching-time identific ation on CR02 CT . The switching-time function t s ( x 0 ) asso ciated with CR02 CT is identified following the pro cedure outlined in Algorithm 1. Although an analytical expression for the switc hing time can b e derived for this example, suc h a closed-form relationship is not generally av ailable in mul- tiparametric dynamic optimization problems. In more complex systems, s witc hing instan ts ma y b e defined im- plicitly through nonlinear junction and jump conditions, making direct analytical deriv ation infeasible. F or this reason, a sampling-based iden tification and fitting strat- egy is employ ed. (1) Identify the switching structure at a represen tative seed. A representativ e initial state x ⋆ 0 ∈ CR02 CT is selected and the con tinuous-time optimalit y condi- tions are solved. The resulting optimal tra jectory exhibits a single switc hing ev en t corresp onding to the release of the low er output constraint. The iden- tified switching structure is therefore S : y min activ e → unconstrained , with one switching even t ( N s = 1). (2) Sample initial states within the domain where the same switching structure is observ ed. A finite set of initial states { x ( r ) 0 } N R r =1 ⊂ CR02 CT 10 T able 3 Explicit con trol and state functions for contin uous-time re- gions. Region Inequalit y Con trol and state functions CR01 CT − 1 . 2707 ≤ x 0 ≤ − 0 . 9323 0 ≤ t ≤ T u ( t ) = ( − x 0 − 1) e t x ( t ) = ( x 0 + 1) e t − 1 CR02 CT − 0 . 9323 ≤ x 0 ≤ − 0 . 5 0 ≤ t ≤ t s u 0 ( t ) = ( − x 0 − 1) e t x 0 ( t ) = ( x 0 + 1) e t − 1 t s ≤ t ≤ T u 1 ( t ) = − 1 2 e − t + t s x 1 ( t ) = − 1 2 e − t + t s CR03 CT − 0 . 5 ≤ x 0 ≤ 0 . 5 0 ≤ t ≤ T u ( t ) = x 0 e − t x ( t ) = x 0 e − t CR04 CT 0 . 5 ≤ x 0 ≤ 0 . 9323 0 ≤ t ≤ t s u 0 ( t ) = (1 − x 0 ) e t x 0 ( t ) = ( x 0 − 1) e t + 1 t s ≤ t ≤ T u 1 ( t ) = 1 2 e − t + t s x 1 ( t ) = 1 2 e − t + t s CR05 CT 0 . 9323 ≤ x 0 ≤ 2 0 ≤ t ≤ T u ( t ) = (1 − x 0 ) e t x ( t ) = ( x 0 − 1) e t + 1 is selected such that the switc hing structure is pre- serv ed for all samples. (2a) Solve the optimality system consistent with the iden tified structure. F or each sample x ( r ) 0 , the dy- namic optimization problem is solved by enforc- ing the same t w o-arc structure as in CR02 CT : an initial arc with active low er output constraint fol- lo w ed by an unconstrained arc. The state, control, adjoin t, and m ultiplier tra jectories are computed accordingly . (2b) Compute the switc hing time for each sample b y en- forcing junction and jump conditions. F or each x ( r ) 0 , the switc hing time t ( r ) s is obtained b y enforcing the junction condition asso ciated with leaving the ac- tiv e constraint. In this example, the condition re- duces to x  t ( r ) s  = − 0 . 5 , whic h yields t ( r ) s = ln  1 2( x ( r ) 0 + 1)  . (2c) Retain only switching times consistent with fea- sibilit y and structural assumptions. The pair  x ( r ) 0 , t ( r ) s  is retained only if the switc hing time satisfies t ( r ) s ∈ [0 , 2] and the constraint activity b efore and after the switching instan t is consistent with the identified structure S . (3) Fit an empirical switc hing-time function. Using the retained sample pairs { ( x ( r ) 0 , t ( r ) s ) } , an empirical ap- T able 4 Switc hing-time functions (exact and fitted) for switching re- gions. Region Exact t s ( x 0 ) Fitted ˆ t s ( x 0 ) R 2 CR02 CT ln  1 2( x 0 +1)  ˆ t s ( x 0 ) = − 25 . 6336 x 3 0 − 46 . 1437 x 2 0 − 30 . 0310 x 0 − 6 . 7059 0.999037 CR04 CT ln  − 1 2( x 0 − 1)  ˆ t s ( x 0 ) = +25 . 6336 x 3 0 − 46 . 1437 x 2 0 + 30 . 0310 x 0 − 6 . 7059 0.999037 pro ximation of the switching-time function is con- structed. A cubic p olynomial is found to provide an accurate representation ov er CR02 CT : ˆ t s ( x 0 ) = − 25 . 63 x 3 0 − 46 . 14 x 2 0 − 30 . 03 x 0 − 6 . 71 , with co efficient of determination R 2 = 0 . 9990. (4) Return the switching-time function asso ciated with the structure. The fitted function ˆ t s ( x 0 ) is asso ci- ated with structure S on CR02 CT . The exact ex- pression is av ailable here for reference: t s ( x 0 ) = ln  1 2( x 0 + 1)  . The resulting exact and fitted switching-time functions obtained from the ab ov e procedure are summarized in T able 4. F or each switc hing region, the exact expression is a v ailable for reference, while the fitted function pro- vides a practical approximation when an explicit ana- lytical form is not av ailable. 4.1.3 Discr ete-Time Counterp art T o compare with the contin uous-time multiparamet- ric solution, the problem is discretized on a uniform grid t k = k h , k = 0 , . . . , N , with step size h = T / N and a zero-order hold (ZOH) input parameterization, u ( t ) = u k , t ∈ [ t k , t k +1 ) , k = 0 , . . . , N − 1 . Under ZOH, the con tinuous-time dynamics ˙ x ( t ) = Ax ( t ) + B u ( t ) admit the exact discrete-time represen- tation x k +1 = A d x k + B d u k , (70) where A d = e Ah and B d = R h 0 e Aτ B dτ , which is the standard mapping used in MPC for linear systems [27]. F or the sp ecific system ˙ x ( t ) = − u ( t ), we hav e A = 0 and B = − 1, hence x k +1 = x k − hu k , k = 0 , . . . , N − 1 , (71) whic h is exact under the ZOH input parameterization. The contin uous-time constraints y ( t ) = x ( t ) + u ( t ) , − 1 ≤ y ( t ) ≤ 1 , u ( t ) ≤ 2 are enforced at the grid no des only , − 1 ≤ x k + u k ≤ 1 , u k ≤ 2 , k = 0 , . . . , N − 1 , (72) 11 T able 5 Represen tative explicit DT critical regions for N = 5 ( h = 0 . 4): affine laws u ( θ ) = K u θ + k u and x ( θ ) = K x θ + k x . Region θ interv al u ( θ ) x ( θ ) CR DT 01 [ − 1 . 52 , − 0 . 89]     − 1 − 1 . 4 − 1 . 96 − 2 . 744 − 3 . 8416     θ +     − 1 − 1 . 4 − 1 . 96 − 2 . 744 − 3 . 8416         1 . 4 1 . 96 2 . 744 3 . 8416 5 . 37824     θ +     0 . 4 0 . 96 1 . 744 2 . 8416 4 . 37824     CR DT 06 [ − 0 . 55 , 0 . 55]     0 . 815 0 . 546 0 . 363 0 . 239 0 . 153     θ     0 . 674 0 . 456 0 . 310 0 . 215 0 . 153     θ CR DT 11 [0 . 89 , 2]     − 1 − 1 . 4 − 1 . 96 − 2 . 744 − 3 . 8416     θ +     1 1 . 4 1 . 96 2 . 744 3 . 8416         1 . 4 1 . 96 2 . 744 3 . 8416 5 . 37824     θ −     0 . 4 0 . 96 1 . 744 2 . 8416 4 . 37824     T able 6 Num b er of DT critical regions versus discretization level N . N Number of DT critical regions 5 11 10 21 15 31 30 61 together with the initial condition x 0 = θ , θ ∈ [ − 2 , 2]. The contin uous-time cost J = 1 2 x ( T ) 2 + 1 2 Z T 0  x ( t ) 2 + u ( t ) 2  dt is approximated by a forward-Euler quadrature, J N = 1 2 x 2 N + 1 2 N − 1 X k =0 h  x 2 k + u 2 k  , (73) whic h yields the discrete-time parametric quadratic pro- gram min { x k ,u k } J N s.t. x k +1 = x k − hu k , k = 0 , . . . , N − 1 , − 1 ≤ x k + u k ≤ 1 , k = 0 , . . . , N − 1 , u k ≤ 2 , k = 0 , . . . , N − 1 , x 0 = θ , θ ∈ [ − 2 , 2] . (74) The explicit multiparametric solution (critical-region partition and affine laws) is computed for (74) using PPOPT . F or N = 5, the parameter space is partitioned in to 11 discrete-time critical regions. T able 5 rep orts three represen tativ e regions from the complete set of 11. Within each discrete-time critical region, the optimal con trol and state tra jectories are given by affine func- tions of θ . The larger num b er of DT regions compared to the contin uous-time partition reflects the gro wth of ad- missible KKT-activ e sets induced b y time discretization and no de-wise constraint enforcement. T o illustrate the impact of discretization on explicit- solution complexit y , T able 6 rep orts the num b er of DT critical regions for increasing N . T able 6 shows that the num b er of DT critical re- gions increases with N , consistent with the combinato- rial growth in admissible active-constrain t patterns as the num b er of grid no des increases. 4.2 Example 2: Two-State System with Mixe d Output and Input Constr aints This example extends the framework to a tw o-state system with path constraints, demonstrating the ap- proac h on a more complex dynamic system. 4.2.1 Pr oblem Statement min u ( · ) 1 2 x ( T ) ⊤ x ( T ) + 1 2 Z T 0  x ( t ) ⊤ x ( t ) + u ( t ) 2  dt s.t. ˙ x ( t ) =  0 − 1 − 1 0  x ( t ) +  1 0  u ( t ) , t ∈ [0 , T ] , y ( t ) =  − 1 1 1 − 1  x ( t ) +  1 − 1  u ( t ) , t ∈ [0 , T ] , y ( t ) ≤  1 . 2 2  , t ∈ [0 , T ] , x (0) = x 0 , − 2 ≤ x 0 ,i ≤ 2 , i = 1 , 2 , T = 2 . (75) The initial state ( x 0 , 1 , x 0 , 2 ) is treated as the parameter in problem (75), and the ob jectiv e is to express the optimal con trol law explicitly as a function of the initial state, sub ject to path constraints enforced ov er the full horizon. 4.2.2 Continuous-Time Multip ar ametric Solution The parameter domain is Θ := { ( x 0 , 1 , x 0 , 2 ) ∈ R 2 | − 2 ≤ x 0 , 1 ≤ 2 , − 2 ≤ x 0 , 2 ≤ 2 } . (1) Initialize the unexplored parameter set: Θ rest ← Θ . (2a) Fix a seed p oint in the initial set: x (2) 0 =  0 0 . 5  ∈ Θ rest . (2b) Solve the dynamic optimization problem at x 0 = x (2) 0 . The optimal structure consists of tw o arcs, Arc A (active) : 0 ≤ t ≤ t s , Arc B (inactive) : t s ≤ t ≤ T , where the first path constraint is active on Arc A, g 1 ( t ) = y 1 ( t ) − 1 . 2 = 0 , µ 1 ( t ) ≥ 0 , and Arc B is unconstrained with µ 1 ( t ) = µ 2 ( t ) = 0. (2c) Closed-form expressions on each arc. Arc B (inactive, t s ≤ t ≤ T ). 12 On this arc, µ 1 ( t ) = µ 2 ( t ) = 0. The stationarity con- dition gives u ( t ) = − λ 1 ( t ). Substituting into the state- costate equations yields the linear homogeneous system d dt    x 1 ( t ) x 2 ( t ) λ 1 ( t ) λ 2 ( t )    = A in    x 1 ( t ) x 2 ( t ) λ 1 ( t ) λ 2 ( t )    , A in =    0 − 1 − 1 0 − 1 0 0 0 − 1 0 0 1 0 − 1 1 0    . A closed-form solution is x 1 ( t ) = √ 2 C 1 e − √ 2 t − √ 2 C 2 e √ 2 t + C 3 e − t + C 4 e t , x 2 ( t ) = C 1 e − √ 2 t + C 2 e √ 2 t + C 3 e − t − C 4 e t , λ 1 ( t ) = C 1 e − √ 2 t + C 2 e √ 2 t , λ 2 ( t ) = C 3 e − t + C 4 e t , t ∈ [ t s , T ] . where C 1 , . . . , C 4 are arc constants. Arc A (constraint 1 active, 0 ≤ t ≤ t s ). On this arc, g 1 ( t ) = 0 and µ 1 ( t ) ≥ 0. Using stationar- it y and differentiating the active constrain t leads to an augmen ted linear system in ( x, λ, µ 1 ), d dt      x 1 ( t ) x 2 ( t ) λ 1 ( t ) λ 2 ( t ) µ 1 ( t )      = A ac      x 1 ( t ) x 2 ( t ) λ 1 ( t ) λ 2 ( t ) µ 1 ( t )      , A ac =      0 − 1 − 1 0 − 1 − 1 0 0 0 0 − 1 0 0 1 1 0 − 1 1 0 − 1 0 1 1 − 1 0      . A closed-form solution is x 1 ( t ) = C ′ 1 e − t − 4 3 C ′ 2 e 2 t , x 2 ( t ) = C ′ 1 e − t + 2 3 C ′ 2 e 2 t − C ′ 3 , λ 1 ( t ) = C ′ 2 e 2 t − C ′ 4 e − 2 t − C ′ 5 e t , λ 2 ( t ) = C ′ 1 e − t − 1 3 C ′ 2 e 2 t − C ′ 3 + C ′ 4 e − 2 t − 2 C ′ 5 e t , µ 1 ( t ) = C ′ 2 e 2 t + C ′ 3 + C ′ 4 e − 2 t + C ′ 5 e t , t ∈ [0 , t s ] . where C ′ 1 , . . . , C ′ 5 are arc constants. T o determine the co efficients in these closed-form ex- pressions, the b oundary conditions at t = 0 and t = T , together with the jump and junction conditions at t = t s , are imp osed on the tw o-arc solution. (i) Boundary conditions: x 1 (0) = x 0 , 1 , x 2 (0) = x 0 , 2 , λ 1 ( T ) = x 1 ( T ) , λ 2 ( T ) = x 2 ( T ) . (ii) Jump (matching) conditions at the switc hing time: x ( t − s ) = x ( t + s ) , λ ( t − s ) = λ ( t + s ) . (iii) Junction condition at the switching time: µ 1 ( t + s ) = 0 , µ 1 ( t − s ) ≥ 0 . T able 7 Selected fitted co efficients used for rep orting in CR02 CT . Co efficien t Fit (in terms of x 0 , 1 , x 0 , 2 , t s ) R 2 C 4 − 0 . 033 x 0 , 1 e − 0 . 88 t s + 0 . 038 x 0 , 2 e − 0 . 73 t s − 0 . 013 0.997 C ′ 4 − 1 . 73 x 0 , 1 e 1 . 61 t s + 1 . 72 x 0 , 2 e 1 . 62 t s 0.999 C ′ 5 0 . 022 x 0 , 1 e − 0 . 88 t s − 0 . 025 x 0 , 2 e − 0 . 73 t s + 0 . 009 0.997 C 3 − 2 . 12 x 0 , 1 e 0 . 87 t s + 3 . 07 x 0 , 2 e 0 . 68 t s 0.999 T able 8 Fitted switching-time functions for the CT switc hing regions in Example 2. Region Fitted ˆ t s ( x 0 , 1 , x 0 , 2 ) R 2 CR02 CT ˆ t s = − 1 . 10 − 4 . 49 x 0 , 1 + 4 . 48 x 0 , 2 − 4 . 68 x 2 0 , 1 + 9 . 35 x 0 , 1 x 0 , 2 − 4 . 67 x 2 0 , 2 − 2 . 17 x 3 0 , 1 + 6 . 49 x 2 0 , 1 x 0 , 2 − 6 . 49 x 0 , 1 x 2 0 , 2 + 2 . 16 x 3 0 , 2 0.999 CR04 CT ˆ t s = − 1 . 13 + 2 . 79 x 0 , 1 − 2 . 79 x 0 , 2 − 1 . 81 x 2 0 , 1 + 3 . 63 x 0 , 1 x 0 , 2 − 1 . 82 x 2 0 , 2 + 0 . 52 x 3 0 , 1 − 1 . 57 x 2 0 , 1 x 0 , 2 + 1 . 57 x 0 , 1 x 2 0 , 2 − 0 . 52 x 3 0 , 2 0.999 After applying these conditions and simplifying, the follo wing explicit relations are obtained: C ′ 3 = − 0 . 6 , C ′ 1 = x 0 , 1 + 2 x 0 , 2 + 2 C ′ 3 3 , C ′ 2 = x 0 , 2 − x 0 , 1 + C ′ 3 2 , and the inactive-arc constants C 1 and C 2 can be written as linear functions of ( C 3 , C 4 ), C 1 = 169 . 21 C 4 − 0 . 81 C 3 , C 2 = 0 . 28 C 4 + 0 . 0028 C 3 . F or the remaining co efficients, the explicit expressions b ecome length y . In general, these coefficients depend on ( x 0 , 1 , x 0 , 2 ) and include exp onential terms inv olving the switc hing time t s ( x 0 , 1 , x 0 , 2 ). F or compact rep orting, we appro ximate the longest expressions with fitted forms that preserv e this dep endence. T able 7 lists the fitted co efficien ts used in CR02 CT . Moreo v er, the switc hing time t s is determined as part of the algebraic solution asso ciated with the tw o-arc structure. F or parameter v alues where switc hing o ccurs, the function t s ( x 0 , 1 , x 0 , 2 ) is obtained by sampling the parameter space and solving the dynamic optimization problem. The resulting fitted switching-time functions are rep orted in T able 8. The resulting con tin uous-time critical regions are summarized in T able 9, together with their corresp ond- ing arc structures and defining inequalities in the pa- rameter space. A graphical represen tation of the resulting critical- region partition is shown in Figure 3. Eac h colored region in Figure 3 corresp onds to one of the active-arc structures listed in T able 9. 4.2.3 Discr ete-Time Counterp art The discrete-time counterpart of Example 2 is ob- tained using the same mo deling and discretization c hoices adopted in Example 1. In the follo wing, w e rep ort discrete-time critical-region maps for tw o dis- 13 T able 9 Activ e-arc structures and critical-region descriptions for Ex- ample 2. Region Arc structure Region description (b orders) CR01 CT unconstrained {− 3 . 36 x 0 , 1 + 3 . 36 x 0 , 2 − 1 . 2 ≤ 0 , 3 . 36 x 0 , 1 − 3 . 36 x 0 , 2 − 2 ≤ 0 } CR02 CT y 1 active → unconstrained {− 3 . 36 x 0 , 1 + 3 . 36 x 0 , 2 − 1 . 2 ≥ 0 , − 1 . 00 x 0 , 1 + x 0 , 2 − 0 . 61 ≤ 0 } CR03 CT y 1 active (full) {− 1 . 00 x 0 , 1 + x 0 , 2 − 0 . 61 ≥ 0 } CR04 CT y 2 active → unconstrained { 3 . 36 x 0 , 1 − 3 . 36 x 0 , 2 − 2 ≥ 0 , 1 . 00 x 0 , 1 − x 0 , 2 − 1 . 01 ≤ 0 } CR05 CT y 2 active (full) { 1 . 00 x 0 , 1 − x 0 , 2 − 1 . 01 ≥ 0 } Fig. 3. Contin uous-time critical-region map for Example 2. T able 10 Represen tative discrete-time critical region for Example 2 (CR01 DT , N = 5). V ariable Index Affine law u ( θ ) k = 0 u 0 = 1 . 00 x 0 , 1 − 1 . 00 x 0 , 2 + 1 . 20 k = 1 u 1 = 1 . 98 x 0 , 1 − 1 . 98 x 0 , 2 + 1 . 79 k = 2 u 2 = 3 . 93 x 0 , 1 − 3 . 93 x 0 , 2 + 2 . 96 k = 3 u 3 = 7 . 81 x 0 , 1 − 7 . 81 x 0 , 2 + 5 . 28 k = 4 u 4 = 15 . 48 x 0 , 1 − 15 . 48 x 0 , 2 + 9 . 89 x ( θ ) k = 1 x 1 , 1 = 1 . 49 x 0 , 1 − 0 . 82 x 0 , 2 + 0 . 49 x 1 , 2 = − 0 . 49 x 0 , 1 + 1 . 16 x 0 , 2 − 0 . 10 k = 2 x 2 , 1 = 2 . 63 x 0 , 1 − 2 . 18 x 0 , 2 + 1 . 31 x 2 , 2 = − 1 . 31 x 0 , 1 + 1 . 75 x 0 , 2 − 0 . 45 k = 3 x 3 , 1 = 5 . 00 x 0 , 1 − 4 . 69 x 0 , 2 + 2 . 82 x 3 , 2 = − 2 . 81 x 0 , 1 + 3 . 11 x 0 , 2 − 1 . 27 k = 4 x 4 , 1 = 9 . 76 x 0 , 1 − 9 . 56 x 0 , 2 + 5 . 74 x 4 , 2 = − 5 . 72 x 0 , 1 + 5 . 92 x 0 , 2 − 2 . 95 k = 5 x 5 , 1 = 19 . 26 x 0 , 1 − 19 . 13 x 0 , 2 + 11 . 48 x 5 , 2 = − 11 . 45 x 0 , 1 + 11 . 59 x 0 , 2 − 6 . 35 cretization levels and pro vide representativ e affine la ws for one region. T able 10 rep orts the affine state and control laws for a representativ e region, CR01 DT , obtained for N = 5. Finally , T able 11 summarizes the num b er of discrete- time critical regions obtained for increasing discretiza- tion levels. (a) N = 5 (b) N = 10 Fig. 4. Discrete-time critical-region maps for Example 2: (a) N = 5 (11 critical regions); (b) N = 10 (21 critical regions). T able 11 Num b er of DT critical regions v ersus discretization level N for Example 2. N Num b er of DT critical regions 5 11 8 17 12 25 20 41 4.3 Comp arison and R emarks Figure 5 reveals the essen tial distinction b etw een the tw o approaches through a representativ e tra jectory from Example 2. The contin uous-time solution exhibits a smo oth exp onential profile with a switching time at t s = 0 . 1396. With only N = 5 interv als, the discrete- time approximation pro duces a crude staircase tra jec- tory that misses b oth the switching instan t and the exp onen tial structure. Refining to N = 20 impro v es ac- curacy considerably , but the n um b er of critical regions gro ws from 11 to 41, whereas the contin uous-time parti- tion remains fixed at 5 regions. This reveals a fundamen- tal trade-off: discrete-time formulations m ust balance appro ximation accuracy against partition complexity . 14 Fig. 5. Optimal con trol tra jectories for Example 2 at x 0 = ( − 0 . 95 , − 1 . 65): contin uous-time (blue) versus discrete– time with N = 5 (orange) and N = 20 (green). This visual difference stems from how each approac h handles control parameterization and time. Discrete- time formulations adopt a zero-order hold assumption, restricting the control to be piecewise constant across N in terv als. This introduces N state and control v ariables as decision v ariables, leading to Karush–Kuhn–T uc ker conditions that form a finite-dimensional algebraic system scaling with N . Contin uous-time form ulations imp ose no such restriction—time remains a contin u- ous v ariable, and the optimal control structure emerges from Pon tryagin’s Maximum Principle through cou- pled state–adjoint differential equations. The problem dimension stays fixed, determined only by the system order. The resulting solutions take distinctly different forms. F or linear–quadratic problems, contin uous-time solu- tions are closed-form exp onen tial functions of time, ex- plicitly dep ending on x 0 , parameter-dep enden t switc h- ing times t s ( x 0 ), and con tinuous time t (T ables 3 and Section 4.2.2). Discrete-time solutions are piece- wise affine in θ at each grid p oint k , taking the form u ( θ ) = K u θ + k u and x ( θ ) = K x θ + k x (T ables 5 and 10), with time entering implicitly through the index rather than explicitly . A particularly v aluable feature of the contin uous-time approac h is its explicit identification of switching times t s ( x 0 ), revealing precisely when the activ e constraint set c hanges as a function of initial conditions. T able 4 pro- vides analytical expressions characterizing these tran- sitions. Discrete-time formulations cannot capture this structure directly—transitions o ccur only at predefined grid p oin ts, making the apparent switching b ehavior de- p enden t on discretization rather than system dynamics. The partition complexity difference is substantial. F or Example 1, the discrete-time form ulation requires 11 re- gions at N = 5 and 61 regions at N = 30 (T ables 6 and 11), while the con tinuous-time partition contains just 5 regions regardless of discretization. This com bina- torial growth arises from the increasing n um b er of con- strain t activ ation patterns across more time nodes, forc- ing practitioners to balance accuracy against computa- tional burden. Constrain t enforcemen t rev eals another key differ- ence. Con tinuous-time formulations enforce path con- strain ts ov er the entire horizon, for all t ∈ [0 , T ], yielding a fixed feasible region. Discrete-time form ulations chec k constrain ts only at grid p oints, p oten tially missing inter- sample violations. In Example 1, the con tinuous-time feasible set is θ ∈ [ − 1 . 27 , 2 . 0], whereas discrete-time with N = 5 admits θ ∈ [ − 1 . 52 , 2 . 0]. As N increases, the discrete-time feasible region con tracts tow ard its con tin uous-time c oun terpart. Despite these differences, both formulations partition the parameter space using the same principle: within eac h critical region, the activ e constraint set remains in v arian t and the corresp onding Lagrange multipliers satisfy strict complemen tarit y . The distinction lies in represen tation—con tinuous-time solutions encode this through switc hing times and exp onential tra jectories, while discrete-time solutions appro ximate it through an expanding collection of affine segments. 5 Discussion The comparison in Section 4 shows that contin uous- time form ulations ac hiev e an order-of-magnitude reduc- tion in critical regions—5 versus 61 in Example 1. This directly impacts online implementation: lo cating the curren t parameter among fewer regions is substantially faster, determining practical viability for systems re- quiring rapid control up dates or op erating under strict computational constraints. The examples focus on linear-quadratic problems, whic h arise in pro cess control, tra jectory planning, and resource allo cation. F or this class, P ontry agin’s Maxi- m um Principle yields analytical solutions. Several direc- tions warran t inv estigation for extending contin uous- time multiparametric metho ds b eyond this setting. When uncertaint y such as disturbances affects the sys- tem, robust contin uous-time multiparametric formula- tions b ecome necessary to maintain feasibility and con- strain t satisfaction. Developing systematic approaches for handling parametric uncertaint y in the contin uous- time framew ork represen ts an imp ortant researc h direc- tion. Computational to ols sp ecifically designed for solving linear-quadratic contin uous-time multiparametric opti- mization problems do not y et exist. Developing such soft- w are that automates the workflo w from problem sp eci- fication to p arametric solution construction would facil- itate practical implementation and wider adoption. Problems with multiple switching times in tro duce sig- nifican t complexity . When the n um b er of switc hes v aries with parameters, critical regions may b ecome noncon- v ex. Additionally , when the minimum o ver adjoin t mul- tipliers or maximum ov er path constraint violations do not o ccur at time endp oints (0 or T ), the critical region b oundaries become nonlinear rather than linear. Sys- tematic methods for handling these cases while preserv- ing computational adv antages merit inv estigation. Extending to nonlinear systems presen ts b oth c hal- lenges and opportunities. Recen t results [1] demonstrate that strict conv exity of the Hamiltonian in the control v ariable ensures global optimality for contin uous-time nonlinear optimal control problems, even when dy- namics and costs are nonconv ex. This con trasts with discrete-time formulations, where Karush-Kuhn-T uck er 15 Fig. 6. Next generation of parametric optimization and con- trol (P AROC) framework. The highlighted path indicates the tra jectory tow ard multiparametric methods in ODE spaces. conditions pro vide only necessary conditions and m ulti- ple local optima may exist. F or nonlinear systems where the Hamiltonian is conv ex in control, Pon tryagin’s Max- im um Principle yields globally optimal solutions—an adv an tage ov er discrete-time approaches. How ever, the resulting b oundary v alue problems rarely admit closed- form solutions, necessitating hybrid numerical-explicit approac hes. Identifying which nonlinear classes p ermit tractable con tinuous-time multiparametric solutions while preserving global optimality guarantees represents a promising research direction. Within the P AROC framework [23] (Figure 6), the con tin uous-time multiparametric approach developed here represents the next generation: working directly with ordinary differen tial equations through P ontry a- gin’s Maxim um Principle rather than algebraic Karush- Kuhn-T uc k er conditions from discretization. Mo del re- duction remains employ ed when needed, but the funda- men tal shift is from algebraic formulations to differen tial equation-based metho ds. F or linear systems, this ODE- based direction is now viable; extending to the broader problem classes outlined ab o ve would establish it as a practical framework for more complex applications. 6 Conclusions This work demonstrates that m ultiparametric op- timal control can b e p erformed directly on ordinary differen tial equations for linear-quadratic systems, ac hieving order-of-magnitude reductions in solution complexit y compared to discrete-time formulations. The contin uous-time partition remains fixed at 5 crit- ical regions regardless of accuracy requirements, while discrete-time partitions grow from 11 to ov er 60 regions as discretization refines—a difference that directly im- pacts b oth offline computation and online region lo okup for real-time mo del predictive control. The computational adv antage stems from w orking directly with P ontry agin’s Maximum Principle rather than discretized Karush-Kuhn-T uck er conditions. By a v oiding artificial time discretization, the con tinuous- time approac h captures the essential switching struc- ture more compactly while providing explicit analytical expressions for switching times t s ( x 0 ) and state-control tra jectories una v ailable from piecewise-affine discrete- time laws. These explicit switching structures also accel- erate gradient-based online optimization by providing analytical information ab out solution top ology . The systematic algorithms presented here—for switc hing structure identification, parametric switch- ing time computation, and critical region construc- tion—establish contin uous-time m ultiparametric pro- gramming as a practical computational framework for linear-quadratic systems. Natural extensions include robust formulations under parametric uncertain ty , non- linear systems where Hamiltonian con vexit y ensures global optimality despite noncon v ex dynamics, and dedicated computational to ols that automate the work- flo w from problem sp ecification to parametric solution construction. By demonstrating viability for the linear- quadratic case, this work op ens a path tow ard the next generation of parametric metho ds op erating directly in ordinary differential equation spaces. References [1] Abhijeet, Mohamed Nav eed Gul Mohamed, Aa yushman Sharma, and Suman Chakrav orty . Conv exity in optimal control problems. In 2024 IEEE 63r d Confer ence on De cision and Control (CDC) , pages 261–266, 2024. [2] Dirk Augustin and Helmut Maurer. Computational sensitivity analysis for state constrained optimal con trol problems. A nnals of Op erations R ese ar ch , 101(1):75–99, 2001. [3] Alb erto Bemp orad, F rancesco Borrelli, and Manfred Morari. Min-max control of constrained uncertain discrete-time linear systems. In Eur op e an Contr ol Conferenc e (ECC) , pages 1651–1656, 2003. [4] Alb erto Bemp orad, Manfred Morari, Vivek Dua, and Efstratios N Pistikopoulos. The explicit linear quadratic regulator for constrained systems. A utomatic a , 38(1):3–20, 2002. [5] Lorenz T Biegler, Stephen L Campb ell, and V olker Mehrmann. Contr ol and optimization with differ ential- algebr aic c onstraints . SIAM, 2012. [6] F rancesco Borrelli, Mato Baoti´ c, Alb erto Bemp orad, and Manfred Morari. Explicit mp c for discrete hybrid systems. IEEE T r ansactions on Automatic Control , 50(12):1919–1925, 2005. [7] Kathryn Eleda Brenan, Stephen L Campb ell, and Linda Ruth Petzold. Numerical solution of initial-value pr oblems in differ ential-algebr aic equations . SIAM, 1995. [8] V assilis M Charitopoulos, Vivek Dua, and Lazaros G Papageorgiou. Closed loop integration of planning, scheduling and control via exact multi-parametric nonlinear programming. In Computer aided chemic al engine ering , volume 40, pages 1273–1278. Elsevier, 2017. [9] V assilis M Charitopoulos, Lazaros G P apageorgiou, and Vivek Dua. Nonlinear mo del-based pro cess operation under uncertainty using exact parametric programming. Engine ering , 3(2):202–213, 2017. [10] Nuno MC De Oliveira and Lorenz T Biegler. An extension of newton-type algorithms for nonlinear process control. Automatic a , 31(2):281–286, 1995. [11] Moritz Diehl, Hans Georg Bo ck, Holger Diedam, and P-B Wieber. F ast direct multiple sho oting algorithms for optimal robot control. In F ast motions in biome chanics and rob otics: optimization and fe e dback c ontr ol , pages 65–93. Springer, 2006. [12] Vladimir Dom brovskii, T atiana Ob yedk o, and Maria Samorodov a. Model predictive control of constrained marko vian jump nonlinear sto chastic systems and p ortfolio 16 optimization under market frictions. Automatic a , 87:61–68, 2018. [13] Paolo F alcone, F rancesco Borrelli, Jahan Asgari, Hongtei Eric Tseng, and Dav or Hrov at. Predictive active steering control for autonomous vehicle systems. IEEE T ransactions on c ontr ol systems technolo gy , 15(3):566–580, 2007. [14] Gregory F ran¸ cois, Bala Sriniv asan, and Dominique Bonvin. Use of measurements for enforcing the necessary conditions of optimality in the presence of constraints and uncertaint y . Journal of Pro c ess Contr ol , 15(6):701–712, 2005. [15] Izrail Moiseevitch Gelfand, Richard A Silv erman, et al. Calculus of variations . Courier Corp oration, 2000. [16] Dustin Kenefak e and Efstratios N Pistikopoulos. Pp opt- multiparametric solver for explicit mp c. In Computer Aide d Chemic al Engine ering , v olume 51, pages 1273–1278. Elsevier, 2022. [17] Richard E Kopp. Pon tryagin maximum principle. In Mathematics in Scienc e and Engine ering , volume 5, pages 255–279. Elsevier, 1962. [18] Konstantinos I Kouramas, Christos Panos, Nuno P F a ´ ısca, and Efstratios N Pistikopoulos. An algorithm for robust explicit/multi-parametric mo del predictive con trol. Automatic a , 49(2):381–389, 2013. [19] David G. Luenberger. Optimization by V e ctor Sp ac e Metho ds . John Wiley & Sons, New Y ork, 1997. [20] Kenneth R Muske and James B Rawlings. Mo del predictive control with linear mo dels. AIChE Journal , 39(2):262–287, 1993. [21] Richard Ob erdieck and Efstratios N Pistikopoulos. Explicit hybrid mo del-predictive con trol: The exact solution. Automatic a , 58:152–159, 2015. [22] Efstratios N Pistikopoulos, Nikolaos A Diangelakis, and Richard Oberdieck. Multi-p ar ametric optimization and c ontr ol . John Wiley & Sons, 2020. [23] Efstratios N Pistikopoulos, Nikolaos A Diangelakis, Richard Oberdieck, Maria M P apathanasiou, Ioana Nascu, and Muxin Sun. Paroc—an integrated framework and softw are platform for the optimisation and adv anced model-based control of process systems. Chemical Engineering Scienc e , 136:115– 138, 2015. [24] Efstratios N Pistikopoulos, Viv ek Dua, Nikolaos A Bozinis, Alberto Bemporad, and Manfred Morari. On- line optimization via off-line parametric optimization to ols. Computers & Chemical Engine ering , 26(2):175–185, 2002. [25] Efstratios N. Pistikopoulos, Michael C. Georgiadis, and Vivek Dua, editors. Multi-Par ametric Mo del-Based Contr ol . Wiley- VCH, W einheim, Germany , 2007. [26] Lev Semenovic h P ontry agin, Vladimir Grigorevic h Bolty anskii, Rev az V alerianovich Gamkrelidze, and Evgenii F rolovic h Mishchenk o. The Mathematic al The ory of Optimal Pro c esses . Wiley , New Y ork, 1962. [27] James Blake Rawlings, David Q Mayne, Moritz Diehl, et al. Mo del pr edictive c ontr ol: the ory, c omputation, and design , volume 2. Nob Hill Publishing Madison, WI, 2020. [28] Jacques Ric halet. Industrial applications of mo del based predictive control. Automatic a , 29(5):1251–1274, 1993. [29] V Sakizlis, JD Perkins, and EN Pistikopoulos. Explicit solutions to optimal control problems for constrained contin uous-time linear systems. IEE Pr o c e edings-Contr ol The ory and Applic ations , 152(4):443–452, 2005. [30] V assilis Sakizlis, Nikolaos MP Kakalis, Vivek Dua, Julian D Perkins, and Efstratios N Pistikopoulos. Design of robust mo del-based controllers via parametric programming. Automatic a , 40(2):189–201, 2004. [31] Balasubramaniam Sriniv asan, Dominique Bonvin, Erik Visser, and Sriniv as Palanki. Dynamic optimi zation of batch processes: Ii. role of measurements in handling uncertaint y . Computers & chemical engine ering , 27(1):27–44, 2003. [32] Balasubramaniam Sriniv asan, Sriniv as Palanki, and Dominique Bonvin. Dynamic optimization of batch pro cesses: I. characterization of the nominal solution. Computers & Chemical Engine ering , 27(1):1–26, 2003. [33] Rob ert F Stengel. Optimal c ontr ol and estimation . Courier Corporation, 1994. [34] Victor M Zav ala and Lorenz T Biegler. The advanced- step nmpc controller: Optimality , stability and robustness. Automatic a , 45(1):86–93, 2009. [35] Victor M Zav ala, Carl D Laird, and Lorenz T Biegler. F ast implementations and rigorous models: Can b oth be accommodated in nmp c? International Journal of R obust and Nonline ar Control: IF AC-Affiliate d Journal , 18(8):800– 815, 2008. 17