Time-Certified and Efficient NMPC via Koopman Operator

Time-Certified and Efficien t NMPC via Ko opman Op erator ⋆ Liang W u ∗ Y unhong Che ∗∗ Bo Y ang ∗∗∗ Kangyu Lin ∗∗∗∗ J´ an Drgo ˇ na ∗ ∗ Johns Hopkins University, MD 21218, USA e-mail: { w liang14, jdr gona1 } @jh.e du. ∗∗ Massachusetts Institute of T e chnolo gy, MA 02139, USA e-mail: yunhche@mit.e du. ∗∗∗ Tsinghua University, Beijing 100084, China e-mail: yang-b21@mails.tsinghua.e du.cn. ∗∗∗∗ Kyoto University, Gr aduate Scho ol of Informatics, JPN e-mail: k-lin@sys.i.kyoto-u.ac.jp. Abstract: Certifying and accelerating execution times of nonlinear mo del predictive con trol (NMPC) implemen tations are t w o core requiremen ts. Execution-time certificate guaran tees that the NMPC controller returns a solution b efore the next sampling time , and ac hieving faster w orst-case and a verage execution times further enables its use in a wider set of applications. Ho wev er, NMPC pro duces a nonlinear program (NLP) for which it is c hallenging to derive its execution time certificates. Our previous w orks, (W u and Braatz, 2025a; W u et al., 2025b) pro vide data-indep endent execution time certificates (certified n umber of iterations) for box- constrained quadratic programs (BoxQP). T o apply the time-certified Bo xQP algorithm (W u et al., 2025b) for state-input constrained NMPC, this pap er i) learns a linear model via Ko opman op erator; ii) prop oses a dynamic-relaxation construction approac h yields a structured BoxQP rather than a general QP; iii) exploits the structure of Bo xQP , where the dimension of the linear system solv ed in each iteration is reduced from 5 N ( n u + n x ) to N n u (where n u , n x , N denote the n umber of inputs, states, and length of prediction horizon), yielding substantial speedups (when n x ≫ n u , as in PDE con trol). Keywor ds: Nonlinear mo del predictiv e control, Execution time certificate, Koopman Op erator. 1. INTRODUCTION Mo del predictive control (MPC) is a model-based optimal con trol method widely used in man ufacturing, energy sys- tems, and robotics. A t each sampling instant, MPC solves an online optimization problem defined by a prediction mo del, constraints, and an ob jective. Tw o key requirements for deploying MPC in pro duction are achieving (i) low av erage execution time, and (ii) a w orst-case execution time that remains below the sampling p eriod. While most works focus on improving a verage exe- cution time F erreau et al. (2014); Stellato et al. (2020); W u and Bemp orad (2023b,a), the more critical requiremen t is certifying worst-case execution time, as reflected in the standard assumption that each MPC optimization m ust b e solved before the next sampling instant. The execution time is computed from the w orst-case n umber of floating-p oin t op erations ([flops]) using the appro ximate relation execution time = total [flops] required b y the algorithm a verage [flops] processed p er second [ s ] . ⋆ This research w as supp orted by the Ralph O’Connor Sustainable Energy Institute at Johns Hopkins Univ ersity . where the denominator depends primarily on the em b ed- ded pro cessor technology . Determining the worst-case total [flops] requires kno wing the worst-case n umber of itera- tions of the optimization algorithm used in the NMPC, pro vided that each iteration has a known and uniform [flops]. Certifying the num ber of iterations of an iterative opti- mization algorithm is an op en challenge in online paramet- ric NMPC. Generally , the required num b er of iterations dep ends on the con vergence speed, and the distance b e- t ween the optimal point and the initial point. The latter is hard to b ound in adv ance in parametric NMPC scenarios, where the problem data v aries with the feedback states at each sampling time. F or example, NMPC sc hemes that reform ulate the NLP as a sequence of QPs, via successive linearization or real-time iteration Gros et al. (2020), result in a Hessian matrix that v aries o ver time. Exe cution time c ertific ate problem has recen tly b ecome an activ e researc h topic in MPC field, see (Rich ter et al., 2011; Giselsson, 2012; Patrinos and Bemp orad, 2013; Cimini and Bemp orad, 2017; Arnstr¨ om and Axehill, 2019; Cimini and Bemp orad, 2019; Arnstr¨ om et al., 2020; Arnstr¨ om and Axehill, 2021; Ok a wa and Nonak a, 2021; W u and Braatz, 2025a,b). Among these works, our prior studies (W u and Braatz, 2025a,b) provide iteration bounds that are b oth simple to compute and data-independent. Th us, w e applied them for certifying execution times of input- constrained NMPC problems, see W u et al. (2024a,b), where the Ko opman framew ork and R T I sc heme are used, respectively . But, this is limited to input-constrained NMPC problems b ecause the execution-time-certified al- gorithm in W u and Braatz (2025a,b) is only for b o x- constrained quadratic programs (BoxQP). Later, Ref. (W u et al., 2025a) provided execution-time certification and in- feasibilit y detection for general QPs, though its execution time certificate imp oses a notable loss in practical sp eed. Our latest work (W u et al., 2025b) developed a predictor- corrector interior-point method (IPM) based BoxQP algo- rithm, whic h preserv es execution time certificate and com- parable computation efficiency at the same time. Based on that, this pap er dev elops a time-certified and efficien t NMPC solution via Ko opman op erator. 1.1 Contributions T o certify execution times for state–input constrained NMPC, this pap er makes the follo wing contributions: i) learns a linear high-dimensional mo del via Ko opman op erator; ii) prop oses a dynamic-relaxation construction approac h that puts ℓ 2 norm of the m ulti-step prediction mo del in to the ob jective (rather than handled as equality constrain ts), thereb y yielding a structured BoxQP rather than a general QP; iii) exploits the structure of BoxQP based on the algo- rithm framework (W u et al., 2025b), where the di- mension of the linear system solv ed in each iteration is reduced from 5( n u + n x ) N p to n u N p (where n u , n x , N p denote the n umber of inputs, states, and length of pre- diction hoirzon), yielding substantial speedups (when n x ≫ n u , as in PDE con trol). 2. PROBLEM F ORMULA TIONS This article considers a nonlinear MPC problem (NMPC) for trac king, as shown in (1), Nonlinear MPC: min N − 1 X k =1 ∥ x k +1 − x r ∥ 2 W x + ∥ u k − u r ∥ 2 W u (1a) + ∥ u k − u k − 1 ∥ 2 W ∆ u s.t. x 0 = x ( t ) , u − 1 = 0 , (1b) x k +1 = f ( x k , u k ) , k = 0 , . . . , N − 1 (1c) − 1 ≤ x k +1 ≤ 1 , k = 0 , . . . , N − 1 (1d) − 1 ≤ u k ≤ 1 , k = 0 , . . . , N − 1 (1e) where x ( t ) is the feedback state at the sampling time t , u k ∈ R n u and x k ∈ R n x denote the con trol input and state at the k th time step, the prediction horizon length is N , x r and u r denote the desired tracking reference signal for the state and con trol input, respectively , and W x ≻ 0, W u ≻ 0, and W ∆ u ≻ 0 are weigh ting matrices for the deviation of the state trac king error, con trol input trac king error, and con trol input increment, resp ectiv ely . Assume that W x , W u , and W ∆ u are diagonal. The equalit y constraint (1c) represen ts the nonlinear discrete-time dynamical system. Without loss of generalit y , the state and control input constrain ts, [ x min , x max ] and [ u min , u max ], are scaled to [ − 1 , 1 ]. NMPC (1) is a nonconv ex nonlinear program, whic h can b e infeasible, can b e sensitive to the noise-contaminated x ( t ), can hav e slow solving speed, and can lac k con vergence guaran tees and execution time certificate. A Ko opman framew ork for MPC (Korda and Mezi´ c, 2018a), which transforms NMPC (1) in to a con vex QP problem via data-driv en Ko opman approximations, allo ws the use of execution-time-certified and computationally efficient QP algorithms for real-time applications. 3. PRELIMINAR Y: KOOPMAN TRANSF ORMS NMPC INTO GENERAL QP The Koopman op erator (Ko opman, 1931) pro vides a glob- ally linear representation of nonlinear dynamics. In prac- tice, the infinite-dimensional Ko opman op erator is trun- cated and approximated using data-driv en Extended Dy- namic Mode Decomposition (EDMD) metho ds (Williams et al., 2015, 2016; Korda and Mezi´ c, 2018b). In EDMD sp ecifically , the set of extended observ ables is designed as the “lifted” mapping,  ψ ( x ) u (0)  , where u (0) denotes the first comp onent of the sequence u and ψ ( x ) ≜  ψ 1 ( x ) , · · · , ψ n ψ ( x )  ⊤ ( n ψ ≫ n x ) is chosen from a ba- sis function, e.g., Radial Basis F unctions (RBFs) used in Korda and Mezi´ c (2018a), instead of directly solving for them via optimization. In particular, the approximate Ko opman op erator identification problem is reduced to a least-squares problem, whic h assumes that the sampled data { ( x j , u j ) , ( x + j , u + j ) } (where j denotes the index of data samples and the superscript + denotes the v alue at the next time step) are collected with the up date mapping  x + j u + j  =  f ( x j , u j (0)) S u j  . Then an appro ximation of the Ko opman op erator, A ≜ [ A B ], can b e obtained by solving J ( A, B ) = min A,B N d X j =1 ∥ ψ ( x + j ) − Aψ ( x j ) − B u j (0) ∥ 2 2 . (2) According to Korda and Mezi´ c (2018a), if the designed lifted mapping ψ ( x ) contains the state x after the re- ordering ψ ( x ) ← [ x ⊤ , ψ ( x )] ⊤ , then C = [ I , 0]. The learned linear Ko opman predictor mo del is given as ψ k +1 = Aψ k + B u k , x k +1 = C ψ k +1 , (3) where ψ k ≜ ψ ( x k ) ∈ R n ψ denotes the lifted state space and with ψ 0 = ψ ( x ( t )). If (3) has a high lifted dimension (for a go od approximation), its use in MPC will not increase the dimension of the resulting QP if the high-dimensional observ ables ψ k are eliminated via col( x 1 , x 2 , · · · , x N ) = E ψ  x ( t )  + F col( u 0 , u 1 , · · · , u N − 1 ) , (4) where E ≜      C A C A 2 . . . C A N      , F ≜     C B 0 · · · 0 C AB C B · · · 0 . . . . . . . . . . . . C A N − 1 B C A N − 2 B · · · C B     . (5) Then, b y embedding (4) into the quadratic ob jectiv e (1a) and the state constraint (1d), NMPC (1) can b e reduced to a compact gener al QP with the decision vector z ≜ col( u 0 , · · · , u N − 1 ) ∈ R N × n u as shown in (6), where ¯ W x ≜ blkdiag( W x , · · · , W x ), ¯ W u ≜ blkdiag( W u , · · · , W u ), ¯ R ≜ blkdiag( ¯ W ∆ u , · · · , ¯ W ∆ u , ¯ W N ∆ u ) ¯ W ∆ u = h 2 W ∆ u − W ∆ u − W ∆ u 2 W ∆ u i and ¯ W N ∆ u = h 2 W ∆ u − W ∆ u − W ∆ u W ∆ u i  and h ≜ F ⊤ ¯ W x ( E ψ ( x ( t )) − ¯ x r ) − ¯ W u ¯ u r ( ¯ x r = repmat( x r ) and ¯ u r = repmat( u r )). Note that the computation of the high-dimensional observ able ψ ( x ( t )) is p erformed once and will not be in v olved in the iterations of QP , which minimizes a side effect of the high- dimensional Ko opman op erator. NMPC → Ko opman-QP: min z z ⊤  F ⊤ ¯ W x F + ¯ W u + ¯ R  z + 2 z ⊤ h (6a) s.t. − 1 ≤ E ψ ( x t ) + F z ≤ 1 (6b) − 1 ≤ z ≤ 1 (6c) The resulting Ko opman-QP (6) has N n u decision v ari- ables and 2 N ( n x + n u ) inequalit y constrain ts, indep enden t of n ψ . It ma y b e infeasible due to inevitable mo deling errors in the data-driven Ko opman appro ximation. By softening the state constraints with slack ϵ and adding a large p enalty term ρ ϵ ϵ 2 to the ob jectiv e, feasibilit y is guaran teed. The condensed Koopman-QP is a general QP without exploitable structure, and currently there is no structure-tailored QP solver for Koopman-MPC to enable fast real-time computation. 4. METHODOLOGY (1): KOOPMAN-BO XQP VIA D YNAMICS-RELAXED CONSTRUCTION T o simultaneously handle p oten tial infeasibilit y and ob- tain faster and certified execution time, this article pro- p oses a dynamics-r elaxe d construction that transforms the Ko opman-MPC problem into an alwa ys-feasible paramet- ric BoxQP . The dynamic r elaxation appr o ach does not strictly enforce the Ko opman prediction mo del (4) but instead incorp orates it into the ob jectiv e through a penalty term, whic h results in a strongly con vex Bo xQP (7), where U ≜ col( u 0 , · · · , u N − 1 ), X ≜ col( x 1 , · · · , x N ), and ρ > 0 is a large p enalty parameter reflecting the confidence in the Ko opman mo del’s accuracy . Nonlinear MPC → Ko opman-BoxQP min U,X ( X − ¯ x r ) ⊤ ¯ W x ( X − ¯ x r ) + ( U − ¯ u r ) ⊤ ¯ W u ( U − ¯ u r ) + U ⊤ ¯ RU + ρ ∥ X − E ψ  x ( t )  − F U ∥ 2 2 s.t. − 1 ≤ U ≤ 1 , − 1 ≤ X ≤ 1 (7) F or simplicit y , w e denote the decision vector z ≜ col( U, X ) ∈ R n ( n = N ( n u + n x )), and the prop osed Ko opman-BoxQP (7) is constructed as shown in (8), min z 1 2 z ⊤ H z + z ⊤ h ( x ( t )) s.t. − 1 ≤ z ≤ 1 (8) where H ≜ ρ  F ⊤ F − F ⊤ − F I  +  ¯ W u + ¯ R ¯ W x  ≻ 0 , (9) h ( x ( t )) ≜ ρ  F ⊤ E − E  ψ ( x ( t )) −  ¯ W u ¯ u r ¯ W x ¯ x r  . (10) The dynamics-relaxed Ko opman-BoxQP (7) can be view ed as an alternativ e approach to softening the state con- strain ts, since the error in satisfying the prediction model can b e equiv alen tly interpreted as a relaxation of the state constrain ts. 5. METHODOLOGY (2): TIME-CER TIFIED AND EFFICIENT IPM ALGORITHM F OR BOX QP According to (Bo yd and V andenberghe, 2004, Ch 5), the Karush–Kuhn–T uck er (KKT) condition of Box-QP (8) is the follo wing nonlinear equations, H z + h ( x ( t )) + γ − θ = 0 , (11a) z + ϕ − 1 n = 0 , (11b) z − ψ + 1 n = 0 , (11c) ( γ , θ , ϕ, ψ ) ≥ 0 , (11d) γ ⊙ ϕ = 0 , (11e) θ ⊙ ψ = 0 , (11f ) where γ , θ are the Lagrangian v ariables of the low er and upp er bound, resp ectiv ely , and ϕ, ψ are the slack v ariables of the low er and upper b ound, resp ectiv ely . ⊙ represen ts the Hadamard pro duct, i.e., γ ⊙ ϕ = col( γ 1 ϕ 1 , γ 2 ϕ 2 , · · · , γ n ϕ n ). P ath-following primal–dual IPMs are categorized into tw o t yp es: fe asible and infe asible , distinguished b y whether the initial point satisfies Eqns. (11a)– (11d). F or the comple- men tarity constraints (11e)– (11f), fe asible path-following IPMs require the initial p oint to lie in a narrow neigh- b orhoo d. T o demonstrate this, let us denote the feasible region b y F , i.e., F = { ( z , γ , θ , ϕ, ψ ) : (11a) − (11c) , ( γ , θ , ϕ, ψ ) ≥ 0 } (12) and the set of strictly feasible p oin ts b y F + ≜ { ( z , γ , θ, ϕ, ψ ) : (11a) − (11c) , ( γ , θ , ϕ, ψ ) > 0 } (13) W e also consider the neighborho od N ( β ) ≜  ( z , γ , θ , ϕ, ψ ) ∈ F + :      γ ⊙ ϕ θ ⊙ ψ  − µ 1 2 n     2 ≤ β µ  (14) where the duality measure µ ≜ γ ⊤ ϕ + θ ⊤ ψ 2 n and β ∈ [0 , 1]. F e asible path-follo wing IPMs require the initial p oint: ( z 0 , γ 0 , θ 0 , ϕ 0 , ψ 0 ) ∈ N ( β ) , (15) and computing suc h a p oin t is typically exp ensive for general strictly con vex QPs. 5.1 Cost-fr e e initialization for F e asible IPMs Inspired by our previous work W u and Braatz (2025a), whic h first sho wed that Box-QP admits cost-free initial- ization for feasible IPMs, we prop ose the following initial- ization to ensure ( z 0 , γ 0 , θ 0 , ϕ 0 , ψ 0 ) ∈ N ( β ). R emark 1. F or h ( x ( t )) = 0, the optimal solution of Box- QP (8) is z ∗ = 0. F or h  = 0, first scale the ob jective as min z 1 2 z ⊤ (2 λH ) z + z ⊤ (2 λh ( x ( t ))) whic h do es not affect the optimal solution and can ensure the initial p oin t lies in N ( β ) if λ ← β √ 2 ∥ h ∥ 2 . Then (11a) is replaced b y 2 λH z + 2 λh ( x ( t )) + γ − θ = 0 the initialization strategy for Bo x-QP (8) z 0 = 0 , γ 0 = 1 n − λh ( x ( t )) , θ 0 = 1 n + λh ( x ( t )) , ϕ 0 = 1 n , ψ 0 = 1 n , (16) whic h clearly places this initial p oin t in N ( β ) by its defi- nition in Eqn. (15) (for example,      γ 0 ⊙ ϕ 0 θ 0 ⊙ ψ 0  − µ 1 2 n     2 = β µ , where µ = 1). In particular, this letter chooses β = 1 4 , then λ = 1 4 √ 2 ∥ h ( x ( t )) ∥ 2 . 5.2 Algor ithm and worst-c ase iter ation c omplexity F or simplicity , w e introduce v ≜ col( γ , θ ) ∈ R 2 n , s ≜ col( ϕ, ψ ) ∈ R 2 n . According to Remark 1, we hav e ( z , v , s ) ∈ N ( β ). Then, all the search directions (∆ z , ∆ v , ∆ s ) (for b oth predictor and corrector steps) are obtained as solutions of the following system of linear equations: (2 λH )∆ z + Ω∆ v = 0 (17a) Ω ⊤ ∆ z + ∆ s = 0 (17b) s ⊙ ∆ v + v ⊙ ∆ s = σµ 1 2 n − v ⊙ s (17c) where Ω = [ I , − I ] ∈ R n × 2 n , σ is chosen 0 in predictor steps and 1 in correctors steps, respectively , and µ ≜ v ⊤ s 2 n denotes the dualit y measure. R emark 2. Eqns. (17a) and (17b) imply that ∆ v ⊤ ∆ s = ∆ v ⊤ ( − Ω ⊤ ∆ z ) = ∆ z ⊤ (2 λH )∆ z ≥ 0 , whic h is critical in the following iteration complexity analysis. By letting ∆ γ = σ µ 1 ϕ − γ + γ ϕ ∆ z , ∆ θ = σ µ 1 ψ − θ − θ ψ ∆ z , ∆ ϕ = − ∆ z , ∆ ψ = ∆ z , (18) Eqn. (17) can be reduced into a more compact system of linear equations,  2 λH + diag  γ ϕ  + diag  θ ψ   ∆ z = σµ  1 ϕ − 1 ψ  + γ − θ . (19) The prop osed feasible adaptive-step predictor-corrector IPM algorithm for Bo x-QP (8) is first describ ed in Al- gorithm 1. The conv ergence analysis and w orst-case iteration com- plexit y of Algorithm 1 follo w directly from our earlier w ork W u et al. (2025b). Algorithm 1 Time-certified predictor-corrector IPM for Bo x-QP (8) Input : Given a strictly feasible initial p oin t ( z 0 , v 0 , s 0 ) ∈ N (1 / 4) from Remark 1 and a desired optimal level ϵ . Then the worst-case iteration bound is N max =  log( 2 n ϵ ) − 2 log  1 − 0 . 2348 √ 2 n   . for k = 0 , 1 , 2 , · · · , N max − 1 do 1. if ( v k ) ⊤ s k ≤ ϵ , then break; 2. Compute the predictor direction (∆ z p , ∆ v p , ∆ s p ) b y solving Eqn. (17) with ( z , s, v ) = ( z k , v k , s k ), σ ← 0, and µ ← µ k = ( v k ) ⊤ s k 2 n (in volving Eqns. (19) and (18)); 3. ∆ µ p ← (∆ v p ) ⊤ ∆ s p 2 n ; 4. α k ← min  1 2 , q µ k 8 ∥ ∆ v p ⊙ ∆ s p − ∆ µ p 1 2 n ∥  ; 5. ˆ z k ← z k + α k ∆ z p , ˆ v k ← v k + α k ∆ v p , ˆ s k ← s k + α k ∆ s p ; 6. Compute the corrector direction (∆ z c , ∆ v c , ∆ s c ) b y solving Eqn. (17) with ( z , v, s ) = ( ˆ z k , ˆ v k , ˆ s k ), σ ← 1, and µ ← ˆ µ k = ( ˆ v k ) ⊤ ˆ s k 2 n (in volving Eqns. (19) and (18)); 7. z k +1 ← ˆ z k + ∆ z c , v k +1 ← ˆ v k + ∆ v c , s k +1 ← ˆ s k + ∆ s c ; end Output: z k +1 . L emma 1. (see (W u et al., 2025b, Thm. 1)) Let { ( z k , v k , s k ) } b e generated by Algorithm 1. Then µ k +1 ≤  1 − 0 . 2348 √ 2 n  2 µ k (20) F urthermore, Algorithm 1 requires at most N max =     log( 2 n ϵ ) − 2 log  1 − 0 . 2348 √ 2 n      . (21) Note that in practice, Algorithm 1 exhibits O ( n 0 . 25 ) or O (log n )-order iteration c omplexit y due to the conserv a- tiv eness of our pro of, see W u et al. (2025b). 5.3 Efficient c omputation of Newton systems A t each iteration of Algorithm 1, Steps 2 and 6 of Algo- rithm 1 compute the predictor and corrector directions, resp ectiv ely , using the same coefficient matrix as given by (according to Eqn. (9)): 2 λH + diag  γ k ϕ k + θ k ψ k  =  ¯ H 11 − 2 λρ F ⊤ − 2 λρ F ¯ H 22  ≻ 0 with ¯ H 11 ≜ 2 λρ F ⊤ F + 2 λ ¯ W u + 2 λ ¯ R + diag γ k 1: n 1 ϕ k 1: n 1 + θ k 1: n 1 ψ k 1: n 1 ! , ¯ H 22 ≜ 2 λρI + 2 λ ¯ W x + diag γ k n 1 +1: n ϕ k n 1 +1: n + θ k n 1 +1: n ψ k n 1 +1: n ! , where ¯ H 11 ≻ 0 ∈ R N n u × N n u , ¯ H 22 ≻ 0 ∈ R n x × n x . By assumption the weigh ting matrix W x is diagonal, thus ¯ W x and ¯ H 22 are also diagonal. R emark 3. By exploiting the diagonal structure of ¯ H 22 and the Sch ur Complemen t, the decomposition cost can be reduced to the element-wise division for and the Cholesky factorization ¯ H − 1 22 with the cost O ( N n x ) Chol  ¯ H 11 − ρ 2 F ⊤ ¯ H − 1 22 F  with the cost O (( N n u ) 3 ) . In summary , the transformation from (17) to (19) reduces the problem dimension from 5 N ( n u + n x ) to N ( n u + n x ), and the subsequen t 2 × 2 block reduction of H in (9) further reduces it to N n u . R emark 4. In most PDE control applications, the num b er of discretized spatial states is significan tly larger than the n umber of con trol inputs, i.e., n x ≫ n u . This is wh y our prop osed Algorithm is not only time-certified but also computationally efficien t in large-scale PDE-MPC applications, as demonstrated in Section 6. 6. NUMERICAL EXAMPLES This section applies our prop osed Ko opman-BoxQP to solv e a PDE-MPC problem with state and control input constrain ts, whic h is a large-sized QP problem with 1040 v ariables and 2080 constraints. The PDE plant under consideration is the nonlinear Korteweg-de V ries (KdV) equation that mo dels the propagation of acoustic wa ves in plasma or shallo w water w av es (Miura, 1976) as ∂ y ( t, x ) ∂ t + y ( t, x ) ∂ y ( t, x ) ∂ x + ∂ 3 y ( t, x ) ∂ x 3 = u ( t, x ) (22) where x ∈ [ − π , π ] is the spatial v ariable. Consider the con- trol input u ( t, x ) = P 4 i =1 u i ( t ) v i ( x ), in which the four co- efficien ts { u i ( t ) } are sub ject to the constraint [ − 1 , 1], and v i ( x ) are predetermined spatial profiles given as v i ( x ) = e − 25( x − m i ) 2 , with m 1 = − π / 2, m 2 = − π / 6, m 3 = π / 6, and m 4 = π / 2. The con trol ob jectiv e is to adjust u i ( t ) so that the spatial profile y ( t, x ) tracks the giv en reference signal. The spatial profile y ( t, x ) is uniformly discretized in to 100 spatial nodes. These 100 discretized nodes serve as the system states, eac h constrained within [ − 1 , 1], while the con trol inputs are likewise bounded within [ − 1 , 1]. With a prediction horizon of N = 10, the resulting MPC problem has 1040 v ariables and 2080 constraints. The cost matrices are chosen as W x = I 100 and W u = 0 . 05 I 4 . The state reference x r ∈ R 100 is a sin usoidal signal ov er a 50 s sim ulation, while the input reference is constan t with u r = 0. Data generation and closed- lo op MPC sim ulation are p erformed using a F ourier-based sp ectral method with a split-step scheme for the nonlinear KdV (22), with a sampling time of ∆ t = 0 . 01. The closed-lo op simulation is configured as follows: (i) Data generation: 1000 tra jectories of length 200 are generated from random linear combinations of four spatial profiles with con trol inputs in [ − 1 , 1]. (ii) Ko opman predictor: The lift ψ includes 100 ph ysical states and 200 th in-plate RBFs ( N ψ = 300). The matrices A ∈ R 300 × 300 and B ∈ R 300 × 4 are iden tified via the Moore–Penrose pseudoin v erse, with C = [ I 100 , 0]. (iii) Ko opman-Bo xQP form ulation: Using ( A, B , C ), the m ulti-step mo del is em b edded into (7) with ρ = 10 2 , producing a Box-QP with 1040 v ariables and 2080 constrain ts. Before applying Algorithm 1 to solv e the resulting large- size Ko opman-BoxQP problem (7) in the closed-loop sim u- lation, w e can kno w in adv ance that the maxim um n umber of iterations for n = 1040 , ϵ = 10 − 6 b y Lemma 1: N max =     log( 2 × 1040 10 − 6 ) − 2 log  1 − 0 . 2348 √ 2 × 1040      = 2079 whic h thus offer a w orst-case execution time certificate. In practice, Algorithm 1 requires an av erage of 72 itera- tions and at most 76 iterations across all sampling instan ts during the closed-lo op simulation. The execution times of Algorithm 1 and other state-of-the-art QP solv ers are rep orted in T ab. 1, whic h sho ws that Algorithm 1 is most efficien t. T able 1. Execution time comparison b et ween Algorithm 1 and other state-of-the-art QP solv ers. QP Solv er W orst-case Execution Time [s] 1 Quadprog 0.1127 OSQP 9 . 7 × 10 − 3 SCS 5 . 8 × 10 − 3 Algorithm 1 2 . 163 × 10 − 3 The closed-lo op results in Fig. 1 show that y ( t, x ) accu- rately trac ks the reference while satisfying b oth state and input b ounds in [ − 1 , 1]. 7. CONCLUSION This paper presents a time-certified and efficient NMPC framew ork based on a dynamic-relaxed Ko opman-BoxQP form ulation that transforms NMPC in to a structured Bo x- QP . Building on the time-certified Box-QP solver in W u et al. (2025b), w e in tro duce a structure-exploiting linear- algebra tec hnique that reduces the dimension of the re- sulting linear system, enabling significant sp eedups. F u- ture work will study the stability and robustness of the dynamics-relaxed Ko opman-BoxQP framework and com- pare it against other state-of-the-art optimization solv ers. DECLARA TION OF GENERA TIVE AI AND AI-ASSISTED TECHNOLOGIES IN THE WRITING PR OCESS During the preparation of this work, the authors used ChatGPT to polish paragraphs. After using this to ol/service, the authors review ed and edited the con ten t as needed and tak e full resp onsibilit y for the conten t of the publication. REFERENCES Arnstr¨ om, D. and Axehill, D. (2019). Exact complexit y certification of a standard primal active-set method for quadratic programming. In Pr o c e e dings of the 58th IEEE Confer enc e on De cision and Contr ol , 4317–4324. Arnstr¨ om, D. and Axehill, D. (2021). A Unifying Complex- it y Certification F ramew ork for Activ e-Set Metho ds for Con vex Quadratic Programming. IEEE T r ansactions on Automatic Contr ol , 67(6), 2758–2770. 1 The execution time results are based on MA TLAB’s Quadprog solver, OSQP’s and SCS’s MA TLAB interface ( https://osqp.org/ and https://www.cvxgrp.org/scs/ ), and C-MEX implementation of Algorithm 1, running on a Mac mini with an Apple M4 Chip (10-core CPU and 16 GB RAM). Fig. 1. Closed-lo op simulation of the nonlinear KdV system with the dynamics-relaxed Ko opman-BoxQP controller trac king a time-v arying spatial profile reference. Left: time evolution of the spatial profile y ( t, x ) and the state constrain ts [ − 1 , 1]. Middle: spatial mean of the y ( t, x ) and the state constraints [ − 1 , 1]. Right: the four control inputs and the con trol input constraints [ − 1 , 1]. Arnstr¨ om, D., Bemp orad, A., and Axehill, D. (2020). Complexit y certification of proximal-point metho ds for n umerically stable quadratic programming. IEEE Con- tr ol Systems L etters , 5(4), 1381–1386. Bo yd, S.P . and V andenberghe, L. (2004). Convex Opti- mization . Cambridge Univ ersity Press, U.K. Cimini, G. and Bemp orad, A. (2017). Exact complexity certification of active-set metho ds for quadratic pro- gramming. IEEE T r ansactions on Automatic Contr ol , 62(12), 6094–6109. Cimini, G. and Bemp orad, A. (2019). Complexity and con vergence certification of a block principal piv oting metho d for box-constrained quadratic programs. Auto- matic a , 100, 29–37. F erreau, H.J., Kirc hes, C., P otschk a, A., Bo c k, H.G., and Diehl, M. (2014). qpOASES: A parametric active-set algorithm for quadratic programming. Mathematic al Pr o gr amming Computation , 6(4), 327–363. Giselsson, P . (2012). Execution time certification for gradien t-based optimization in model predictiv e con trol. In Pr o c e e dings of the 51st IEEE Confer enc e on De cision and Contr ol , 3165–3170. Gros, S., Zanon, M., Quirynen, R., Bemp orad, A., and Diehl, M. (2020). F rom linear to nonlinear MPC: Bridg- ing the gap via the real-time iteration. International Journal of Contr ol , 93(1), 62–80. Ko opman, B.O. (1931). Hamiltonian systems and trans- formation in Hilb ert spac e. Pr o c e e dings of the National A c ademy of Scienc es , 17(5), 315–318. Korda, M. and Mezi´ c, I. (2018a). Linear predictors for nonlinear dynamical systems: Ko opman op erator meets mo del predictive con trol. Automatic a , 93, 149–160. Korda, M. and Mezi ´ c, I. (2018b). On conv ergence of extended dynamic mo de decomposition to the Koopman op erator. Journal of Nonline ar Scienc e , 28(2), 687–710. Miura, R.M. (1976). The Korteweg–deVries equation: A surv ey of results. SIAM R eview , 18(3), 412–459. Ok a wa, I. and Nonak a, K. (2021). Linear complementarit y mo del predictiv e control with limited iterations for b ox- constrained problems. A utomatic a , 125, 109429. P atrinos, P . and Bemp orad, A. (2013). An accelerated dual gradient-pro jection algorithm for embedded linear mo del predictive control. IEEE T r ansactions on A uto- matic Contr ol , 59(1), 18–33. Ric hter, S., Jones, C.N., and Morari, M. (2011). Computa- tional complexity certification for real-time MPC with input constraints based on the fast gradien t metho d. IEEE T r ansactions on A utomatic Contr ol , 57(6), 1391– 1403. Stellato, B., Banjac, G., Goulart, P ., Bemporad, A., and Bo yd, S. (2020). OSQP: An operator splitting solver for quadratic programs. Mathematic al Pr o gr amming Computation , 12(4), 637–672. Williams, M.O., Hemati, M.S., Da wson, S.T.M., Kevrekidis, I.G., and Ro wley , C.W. (2016). Extending data-driv en Ko opman analysis to actuated systems. IF A C-Pap ersOnLine , 49(18), 704–709. Williams, M.O., Kevrekidis, I.G., and Rowley , C.W. (2015). A data-driv en approximation of the Ko op- man operator: Extending dynamic mode decomposition. Journal of Nonline ar Scienc e , 25(6), 1307–1346. W u, L. and Bemp orad, A. (2023a). A construction-free co ordinate-descen t augmented-Lagrangian method for em b edded linear MPC based on ARX models. IF AC- Pap ersOnLine , 56(2), 9423–9428. W u, L., Ganko, K., and Braatz, R.D. (2024a). Time- certified Input-constrained NMPC via Koopman op era- tor. IF AC-Pap ersOnLine , 58(18), 335–340. doi:h ttps:// doi.org/10.1016/j.ifacol.2024.09.052. 8th IF A C Confer- ence on Nonlinear Model Predictive Con trol NMPC 2024. W u, L., Xiao, W., and Braatz, R.D. (2025a). EIQP: Execution-time-certified and Infeasibility-detecting QP solv er. IEEE T r ansactions on Automatic Contr ol , 1–16. DOI: 10.1109/TAC.2025.3631575 . W u, L. and Bemp orad, A. (2023b). A Simple and Fast Co ordinate-Descen t Augmented-Lagrangian Solver for Mo del Predictiv e control. IEEE T r ansactions on Auto- matic Contr ol , 68(11), 6860–6866. W u, L. and Braatz, R.D. (2025a). A Direct Optimization Algorithm for Input-Constrained MPC. IEEE T r ansac- tions on A utomatic Contr ol , 70(2), 1366–1373. W u, L. and Braatz, R.D. (2025b). A Quadratic Program- ming Algorithm with O ( n 3 ) Time Complexit y . arXiv pr eprint arXiv:2507.04515 . W u, L., Che, Y., Braatz, R.D., and Drgona, J. (2025b). A Time-certified Predictor-corrector IPM Algorithm for Bo x-QP. IEEE Contr ol Systems L etters , 9, 3059–3064. W u, L., Gank o, K., W ang, S., and Braatz, R.D. (2024b). An Execution-time-certified Riccati-based IPM Algo- rithm for R TI-based Input-constrained NMPC. In 2024 IEEE 63r d Confer enc e on De cision and Contr ol (CDC) , 5539–5545. IEEE.