MPC-RL-MOBO: MPC, RL & Multi-Objective BO

• MPC-RL-MOBO is a framework combining MPC, RL, and MOBO for safe and sample-efficient control in nonlinear, safety-critical systems.
• It parameterizes MPC policies and uses KKT-based implicit differentiation with CDPG for accurate gradient and performance estimation.
• The framework leverages EHVI to balance a four-objective vector, optimizing performance, convergence, critic stability, and Lyapunov-based safety.

Searching arXiv for the target paper and a few adjacent strands to ground the article. {"query":"arXiv (Esfahani et al., 14 Jul 2025) MPC-RL-MOBO differentiable MPC safe Bayesian optimization deterministic policy gradient", "max_results": 10} {"query":"differentiable MPC reinforcement learning arXiv safe Bayesian optimization EHVI", "max_results": 10} MPC-RL-MOBO denotes a framework at the intersection of Model Predictive Control (MPC), Reinforcement Learning (RL), and Multi-Objective Bayesian Optimization (MOBO) for intelligent control of industrial processes. It uses a parameterized MPC problem as the policy class, estimates noisy closed-loop performance and gradient information from rollouts through a Compatible Deterministic Policy Gradient (CDPG) construction, and updates MPC parameters through MOBO with an Expected Hypervolume Improvement (EHVI) acquisition function. The stated purpose is safe, interpretable, and sample-efficient learning under nonlinear dynamics, model mismatch, and safety-critical state and input constraints, while avoiding several limitations attributed to standard MPC-RL approaches, including slow convergence, suboptimal policy learning due to limited parameterization, and safety issues during online adaptation (Esfahani et al., 14 Jul 2025).

1. Problem setting and conceptual scope

Industrial process control in the framework is posed for nonlinear, uncertain plants with safety-critical constraints. The real plant may be represented either as a Markov decision process,

$x_{k+1} \sim P[\cdot \mid x_k, u_k],$

or in control form as

$x_{k+1} = f(x_k, u_k, d_k),$

where $d_k$ denotes disturbances or uncertainties. The nominal model used inside MPC is explicitly imperfect:

$x_{k+1} = \hat f(x_k, u_k) + \delta(x_k, u_k),$

with $\delta$ capturing model mismatch, structural errors, and unmodeled dynamics.

The control problem includes pure input constraints, mixed state-input constraints, and terminal constraints:

$$g(u_k) \le 0, \qquad h(x_k, u_k) \le 0, \qquad h^f(x_{k+N}) \le 0.$$

These may be relaxed through slack variables with large penalties so that feasibility can be maintained even under uncertainty.

Within this setting, the learning problem is to adapt the MPC parameter vector $\theta$ online so as to optimize a multi-objective vector $f(\theta)$ while preserving safety and stabilizing adaptation. The objectives are built from noisy rollout evaluations of the MPC-induced closed-loop behavior together with gradient information obtained from a CDPG critic compatible with the MPC policy. This suggests that the framework treats learning not as direct synthesis of a policy from a generic function approximator, but as structured tuning of an optimizer-defined control law.

2. Parameterized MPC policy and differentiability structure

At each time step, the controller solves a finite-horizon optimal control problem with decision variables $\{\hat x_i,\hat u_i,\eta_i\}_{i=0}^{N}$ and cost

$$J_{\mathrm{MPC}}(x_k;\theta) = \gamma^N \big(V^f_\theta(\hat x_N)+\Gamma_f^\top \eta_N\big) + \sum_{i=0}^{N-1}\gamma^i\big(\ell_\theta(\hat x_i,\hat u_i)+\Gamma^\top \eta_i\big),$$

subject to

$x_{k+1} = f(x_k, u_k, d_k),$0

$x_{k+1} = f(x_k, u_k, d_k),$1

The discount factor satisfies $x_{k+1} = f(x_k, u_k, d_k),$2, and the large penalty weights $x_{k+1} = f(x_k, u_k, d_k),$3 discourage slacks unless they are required for feasibility.

The induced policy is the first optimizer action,

$x_{k+1} = f(x_k, u_k, d_k),$4

The parameter vector $x_{k+1} = f(x_k, u_k, d_k),$5 may include cost weights in $x_{k+1} = f(x_k, u_k, d_k),$6 and $x_{k+1} = f(x_k, u_k, d_k),$7, model parameters or hyperparameters in $x_{k+1} = f(x_k, u_k, d_k),$8, constraint tightenings in $x_{k+1} = f(x_k, u_k, d_k),$9 and $d_k$0, slack penalties $d_k$1, tube or robust design parameters, mixed or terminal constraint parameterizations, and potentially the horizon length $d_k$2, although $d_k$3 is often fixed for differentiability and practical reasons.

A central technical feature is differentiability of $d_k$4 with respect to $d_k$5. Under local uniqueness of the MPC optimum and KKT regularity conditions such as LICQ and SOSC, implicit differentiation of the KKT system yields local sensitivities $d_k$6. The policy is piecewise differentiable across active-set changes, and CDPG together with differentiable MPC uses KKT-based implicit differentiation to compute $d_k$7. In practice, differentiable MPC solvers provide these sensitivities by automatic differentiation and linear algebra on the KKT Jacobian.

This architecture is presented as an interpretable alternative to DNN-based RL methods, with lower computational complexity and greater transparency. The interpretability arises because policy updates operate through explicit MPC quantities—cost weights, constraint parameterizations, nominal model components, and robust design parameters—rather than through opaque latent parameters.

3. RL objective, stage cost, and Compatible Deterministic Policy Gradient

The RL objective is the expected cumulative stage cost generated by the closed-loop policy on the real plant:

$d_k$8

with

$d_k$9

Safety enters the stage cost through an indicator penalty:

$x_{k+1} = \hat f(x_k, u_k) + \delta(x_k, u_k),$0

where

$x_{k+1} = \hat f(x_k, u_k) + \delta(x_k, u_k),$1

Practical evaluations use finite-horizon rollouts of length $x_{k+1} = \hat f(x_k, u_k) + \delta(x_k, u_k),$2, producing noisy samples

$x_{k+1} = \hat f(x_k, u_k) + \delta(x_k, u_k),$3

where $x_{k+1} = \hat f(x_k, u_k) + \delta(x_k, u_k),$4 approximates $x_{k+1} = \hat f(x_k, u_k) + \delta(x_k, u_k),$5 through a finite sum and $x_{k+1} = \hat f(x_k, u_k) + \delta(x_k, u_k),$6 represents evaluation noise from process disturbances, random initial conditions, solver tolerances, and related effects.

For deterministic policies, the framework invokes the deterministic policy gradient theorem:

$x_{k+1} = \hat f(x_k, u_k) + \delta(x_k, u_k),$7

To reduce variance while preserving unbiasedness, it uses a compatible function approximator

$x_{k+1} = \hat f(x_k, u_k) + \delta(x_k, u_k),$8

Its advantage term,

$x_{k+1} = \hat f(x_k, u_k) + \delta(x_k, u_k),$9

targets $\delta$0 and vanishes on-policy at $\delta$1. The baseline $\delta$2 approximates $\delta$3.

The critic parameters $\delta$4 are estimated by Least Squares Temporal Difference methods. The baseline approximation is written as

$\delta$5

where $\delta$6 is the MPC optimal cost-to-go and its sensitivity is available through Lagrangian sensitivity:

$\delta$7

For the action-value approximation, the framework uses LSTDQ on off-policy pairs $\delta$8 with

$\delta$9

and moment condition

$$g(u_k) \le 0, \qquad h(x_k, u_k) \le 0, \qquad h^f(x_{k+N}) \le 0.$$0

Because $$g(u_k) \le 0, \qquad h(x_k, u_k) \le 0, \qquad h^f(x_{k+N}) \le 0.$$1, the batch least-squares solution is

$$g(u_k) \le 0, \qquad h(x_k, u_k) \le 0, \qquad h^f(x_{k+N}) \le 0.$$2

with

$$g(u_k) \le 0, \qquad h(x_k, u_k) \le 0, \qquad h^f(x_{k+N}) \le 0.$$3

$$g(u_k) \le 0, \qquad h(x_k, u_k) \le 0, \qquad h^f(x_{k+N}) \le 0.$$4

For the baseline, LSTDV gives

$$g(u_k) \le 0, \qquad h(x_k, u_k) \le 0, \qquad h^f(x_{k+N}) \le 0.$$5

with

$$g(u_k) \le 0, \qquad h(x_k, u_k) \le 0, \qquad h^f(x_{k+N}) \le 0.$$6

$$g(u_k) \le 0, \qquad h(x_k, u_k) \le 0, \qquad h^f(x_{k+N}) \le 0.$$7

Policy sensitivity is computed by differentiating the MPC KKT system. If $$g(u_k) \le 0, \qquad h(x_k, u_k) \le 0, \qquad h^f(x_{k+N}) \le 0.$$8 is the MPC Lagrangian and $$g(u_k) \le 0, \qquad h(x_k, u_k) \le 0, \qquad h^f(x_{k+N}) \le 0.$$9 collects KKT stationarity and constraint residuals, then under implicit-function-theorem conditions

$\theta$0

Noise in gradient estimates is mitigated through batch averaging, LSTD closed-form solutions, and optional regularization of $\theta$1 and $\theta$2.

4. Multi-objective formulation and Bayesian surrogate modeling

The learning layer optimizes a four-dimensional objective vector

$\theta$3

The objectives are defined as follows.

Objective	Definition	Role stated in the framework
$\theta$4	$\theta$5	finite-horizon approximation of closed-loop performance
$\theta$6	$\theta$7	necessary optimality condition via CDPG gradient magnitude
$\theta$8	$\theta$9 with $f(\theta)$0	critic stability term discouraging growth of the baseline across episodes
$f(\theta)$1	$f(\theta)$2, $f(\theta)$3	Lyapunov-informed stability measure using the MPC value as a surrogate Lyapunov function

Pareto dominance is defined in the standard componentwise minimization form: $f(\theta)$4 dominates $f(\theta)$5 if $f(\theta)$6 for all $f(\theta)$7 and strict inequality holds for at least one objective. The Pareto front is the set of non-dominated points in objective space.

Each objective is modeled by an independent Gaussian process prior

$f(\theta)$8

typically with $f(\theta)$9 and squared exponential kernel

$\{\hat x_i,\hat u_i,\eta_i\}_{i=0}^{N}$0

Observations are noisy,

$\{\hat x_i,\hat u_i,\eta_i\}_{i=0}^{N}$1

and the posterior mean and variance at a test point $\{\hat x_i,\hat u_i,\eta_i\}_{i=0}^{N}$2 are

$\{\hat x_i,\hat u_i,\eta_i\}_{i=0}^{N}$3

$\{\hat x_i,\hat u_i,\eta_i\}_{i=0}^{N}$4

Although derivative information can be incorporated into a joint GP through covariance blocks

$\{\hat x_i,\hat u_i,\eta_i\}_{i=0}^{N}$5

with entries such as

$\{\hat x_i,\hat u_i,\eta_i\}_{i=0}^{N}$6

the framework instead treats gradient information as a separate objective, namely $\{\hat x_i,\hat u_i,\eta_i\}_{i=0}^{N}$7. This is an explicit modeling choice rather than an omission.

5. EHVI acquisition and the role of safety

The acquisition layer is based on Expected Hypervolume Improvement. If $\{\hat x_i,\hat u_i,\eta_i\}_{i=0}^{N}$8 is the current set of non-dominated objective vectors in $\{\hat x_i,\hat u_i,\eta_i\}_{i=0}^{N}$9 with $$J_{\mathrm{MPC}}(x_k;\theta) = \gamma^N \big(V^f_\theta(\hat x_N)+\Gamma_f^\top \eta_N\big) + \sum_{i=0}^{N-1}\gamma^i\big(\ell_\theta(\hat x_i,\hat u_i)+\Gamma^\top \eta_i\big),$$0, and $$J_{\mathrm{MPC}}(x_k;\theta) = \gamma^N \big(V^f_\theta(\hat x_N)+\Gamma_f^\top \eta_N\big) + \sum_{i=0}^{N-1}\gamma^i\big(\ell_\theta(\hat x_i,\hat u_i)+\Gamma^\top \eta_i\big),$$1 is a reference point dominated by all feasible outcomes, then the hypervolume of a set $$J_{\mathrm{MPC}}(x_k;\theta) = \gamma^N \big(V^f_\theta(\hat x_N)+\Gamma_f^\top \eta_N\big) + \sum_{i=0}^{N-1}\gamma^i\big(\ell_\theta(\hat x_i,\hat u_i)+\Gamma^\top \eta_i\big),$$2 relative to $$J_{\mathrm{MPC}}(x_k;\theta) = \gamma^N \big(V^f_\theta(\hat x_N)+\Gamma_f^\top \eta_N\big) + \sum_{i=0}^{N-1}\gamma^i\big(\ell_\theta(\hat x_i,\hat u_i)+\Gamma^\top \eta_i\big),$$3 is $$J_{\mathrm{MPC}}(x_k;\theta) = \gamma^N \big(V^f_\theta(\hat x_N)+\Gamma_f^\top \eta_N\big) + \sum_{i=0}^{N-1}\gamma^i\big(\ell_\theta(\hat x_i,\hat u_i)+\Gamma^\top \eta_i\big),$$4. Sampling $$J_{\mathrm{MPC}}(x_k;\theta) = \gamma^N \big(V^f_\theta(\hat x_N)+\Gamma_f^\top \eta_N\big) + \sum_{i=0}^{N-1}\gamma^i\big(\ell_\theta(\hat x_i,\hat u_i)+\Gamma^\top \eta_i\big),$$5 yields the improvement

$$J_{\mathrm{MPC}}(x_k;\theta) = \gamma^N \big(V^f_\theta(\hat x_N)+\Gamma_f^\top \eta_N\big) + \sum_{i=0}^{N-1}\gamma^i\big(\ell_\theta(\hat x_i,\hat u_i)+\Gamma^\top \eta_i\big),$$6

and EHVI is

$$J_{\mathrm{MPC}}(x_k;\theta) = \gamma^N \big(V^f_\theta(\hat x_N)+\Gamma_f^\top \eta_N\big) + \sum_{i=0}^{N-1}\gamma^i\big(\ell_\theta(\hat x_i,\hat u_i)+\Gamma^\top \eta_i\big),$$7

For independent Gaussian posteriors, EHVI admits closed-form or factorized approximations for small $$J_{\mathrm{MPC}}(x_k;\theta) = \gamma^N \big(V^f_\theta(\hat x_N)+\Gamma_f^\top \eta_N\big) + \sum_{i=0}^{N-1}\gamma^i\big(\ell_\theta(\hat x_i,\hat u_i)+\Gamma^\top \eta_i\big),$$8, and efficient algorithms decompose the integration into orthants defined by the current Pareto set and evaluate expected improvements using Gaussian CDFs. Operationally, EHVI favors points with both attractive posterior means and sufficiently large posterior uncertainty, thereby balancing exploitation and exploration across objectives.

A common misconception is that the BO layer alone is responsible for safety. In the framework, safety is enforced at two layers.

First, the MPC itself contains mixed and terminal constraints together with $$J_{\mathrm{MPC}}(x_k;\theta) = \gamma^N \big(V^f_\theta(\hat x_N)+\Gamma_f^\top \eta_N\big) + \sum_{i=0}^{N-1}\gamma^i\big(\ell_\theta(\hat x_i,\hat u_i)+\Gamma^\top \eta_i\big),$$9 slack variables and large penalties $x_{k+1} = f(x_k, u_k, d_k),$00, which practically enforce constraint satisfaction and feasibility under disturbances. Robust or tube MPC variants can be embedded by parameterizing constraint tightenings and tube sets in $x_{k+1} = f(x_k, u_k, d_k),$01. Stability is further promoted through the objective $x_{k+1} = f(x_k, u_k, d_k),$02, which penalizes positive increments of the MPC value along trajectories and encourages a descent condition of the form

$x_{k+1} = f(x_k, u_k, d_k),$03

Second, Safe BO is described as an optional GP-based extension. If outcome constraints are themselves modeled by GPs, a conservative feasible set may be defined as

$x_{k+1} = f(x_k, u_k, d_k),$04

where $x_{k+1} = f(x_k, u_k, d_k),$05 models rollout-level constraint violations and $x_{k+1} = f(x_k, u_k, d_k),$06 are GP posterior statistics. Acquisition optimization is then restricted to $x_{k+1} = f(x_k, u_k, d_k),$07, with SafeOpt-style expansion only where confidence indicates safety. The paper explicitly states that, in the implementation, safety emphasis is primarily delivered through MPC penalties and the stability objective, while GP-based explicit safety constraints are a natural extension (Esfahani et al., 14 Jul 2025).

An informal Lyapunov proof sketch is given under assumptions that $x_{k+1} = f(x_k, u_k, d_k),$08 is positive definite around a target set and proper, that the closed-loop dynamics under $x_{k+1} = f(x_k, u_k, d_k),$09 are Lipschitz, and that parameter updates keep $x_{k+1} = f(x_k, u_k, d_k),$10 sufficiently small so that

$x_{k+1} = f(x_k, u_k, d_k),$11

for some class-$x_{k+1} = f(x_k, u_k, d_k),$12 function $x_{k+1} = f(x_k, u_k, d_k),$13. The claimed implication is practical stability rather than a full formal theorem.

6. Full algorithmic workflow and computational characteristics

The algorithm begins from an initial safe $x_{k+1} = f(x_k, u_k, d_k),$14, such as conservative cost weights and robust tightenings that make the MPC feasible and stable. It initializes independent GPs for $x_{k+1} = f(x_k, u_k, d_k),$15, sets a reference point $x_{k+1} = f(x_k, u_k, d_k),$16 for EHVI, and chooses the rollout horizon $x_{k+1} = f(x_k, u_k, d_k),$17, discount $x_{k+1} = f(x_k, u_k, d_k),$18, penalty weights $x_{k+1} = f(x_k, u_k, d_k),$19, and the coefficient $x_{k+1} = f(x_k, u_k, d_k),$20 used in $x_{k+1} = f(x_k, u_k, d_k),$21.

Each learning episode $x_{k+1} = f(x_k, u_k, d_k),$22 then follows the sequence below:

• Rollout: apply the MPC policy $x_{k+1} = f(x_k, u_k, d_k),$23 on the real plant for $x_{k+1} = f(x_k, u_k, d_k),$24 steps and record states and inputs.
• Critic updates: evaluate $x_{k+1} = f(x_k, u_k, d_k),$25 and $x_{k+1} = f(x_k, u_k, d_k),$26 from MPC solves and Lagrangian sensitivities; solve LSTDV for $x_{k+1} = f(x_k, u_k, d_k),$27 and LSTDQ for $x_{k+1} = f(x_k, u_k, d_k),$28.
• Gradient estimate: compute $x_{k+1} = f(x_k, u_k, d_k),$29 via KKT-based implicit differentiation; form

$x_{k+1} = f(x_k, u_k, d_k),$30

and average to obtain $x_{k+1} = f(x_k, u_k, d_k),$31.

• Objective evaluations: compute $x_{k+1} = f(x_k, u_k, d_k),$32 from the finite discounted sum of $x_{k+1} = f(x_k, u_k, d_k),$33, $x_{k+1} = f(x_k, u_k, d_k),$34 from episode-level growth of $x_{k+1} = f(x_k, u_k, d_k),$35, and $x_{k+1} = f(x_k, u_k, d_k),$36 from the Lyapunov penalty on $x_{k+1} = f(x_k, u_k, d_k),$37.
• GP update: augment the training set for each objective with $x_{k+1} = f(x_k, u_k, d_k),$38, refit hyperparameters if needed by maximizing marginal likelihood, and update posterior means and variances.
• Acquisition optimization: solve

$x_{k+1} = f(x_k, u_k, d_k),$39

subject to bounds and safety restrictions.

• Termination: stop when EHVI falls below a threshold, the Pareto front stabilizes, or the episode budget is exhausted.

The framework’s sample-efficiency claims are tied to three ingredients: an MPC policy with a small number of interpretable parameters, low-variance gradient estimation via CDPG and LSTD, and MOBO with EHVI for global guidance across multiple objectives. Formal convergence guarantees are not derived. Under standard BO assumptions, EHVI-based MOBO is described as converging to the Pareto set in the limit, with hypervolume improving monotonically in expectation, while differentiable MPC sensitivity and policy-gradient convergence remain local around regular points.

Per-iteration computational costs are stated componentwise. A dense QP or NLP solve with $x_{k+1} = f(x_k, u_k, d_k),$40 decision variables and $x_{k+1} = f(x_k, u_k, d_k),$41 constraints has worst-case cost $x_{k+1} = f(x_k, u_k, d_k),$42, with horizon $x_{k+1} = f(x_k, u_k, d_k),$43 increasing problem size. Sensitivity evaluation requires one linear solve on the KKT system with similar worst-case cubic scaling. LSTD updates cost $x_{k+1} = f(x_k, u_k, d_k),$44 to form $x_{k+1} = f(x_k, u_k, d_k),$45 and $x_{k+1} = f(x_k, u_k, d_k),$46 for inversion. GP updates cost $x_{k+1} = f(x_k, u_k, d_k),$47 per objective with $x_{k+1} = f(x_k, u_k, d_k),$48 training points. EHVI maximization is a nonconvex global optimization over $x_{k+1} = f(x_k, u_k, d_k),$49 and typically uses multi-start plus local gradient-based refinement. Memory is dominated by GP training data and any reused MPC solver artifacts such as KKT factors.

7. Nonlinear CSTR example, practical guidance, and stated limitations

The numerical example uses a nonlinear Continuous Stirred Tank Reactor with reaction $x_{k+1} = f(x_k, u_k, d_k),$50, states

$x_{k+1} = f(x_k, u_k, d_k),$51

and controls

$x_{k+1} = f(x_k, u_k, d_k),$52

The continuous-time dynamics are

$x_{k+1} = f(x_k, u_k, d_k),$53

$x_{k+1} = f(x_k, u_k, d_k),$54

$x_{k+1} = f(x_k, u_k, d_k),$55

The parameters are

$x_{k+1} = f(x_k, u_k, d_k),$56

$x_{k+1} = f(x_k, u_k, d_k),$57

$x_{k+1} = f(x_k, u_k, d_k),$58

with

$x_{k+1} = f(x_k, u_k, d_k),$59

To induce mismatch, the MPC model uses coefficients scaled by factors $x_{k+1} = f(x_k, u_k, d_k),$60. Constraints are

$x_{k+1} = f(x_k, u_k, d_k),$61

$x_{k+1} = f(x_k, u_k, d_k),$62

The MPC configuration uses horizon $x_{k+1} = f(x_k, u_k, d_k),$63, sampling time $x_{k+1} = f(x_k, u_k, d_k),$64 min with RK4 discretization, and exact slack penalties

$x_{k+1} = f(x_k, u_k, d_k),$65

The reference values are

$x_{k+1} = f(x_k, u_k, d_k),$66

$x_{k+1} = f(x_k, u_k, d_k),$67

The RL stage cost is

$x_{k+1} = f(x_k, u_k, d_k),$68

with additional constraint penalties $x_{k+1} = f(x_k, u_k, d_k),$69.

Learning proceeds for 600 episodes, each of duration 3 min or 60 steps, from random initial conditions. The reported objectives are exactly the four quantities $x_{k+1} = f(x_k, u_k, d_k),$70 through $x_{k+1} = f(x_k, u_k, d_k),$71 introduced above. The results reported for this example are threefold. First, MPC-RL-MOBO rapidly expands the Pareto hypervolume and reaches near the optimal performance achieved by a perfect-model MPC; single-objective BO converges slightly slower and less accurately; MPC-RL alone requires approximately 350 episodes to approach the optimum with lower accuracy. Second, the advantage function $x_{k+1} = f(x_k, u_k, d_k),$72 converges toward zero faster under MOBO, which is interpreted as improved satisfaction of necessary optimality conditions. Third, state and input trajectories show that learning tunes MPC parameters from imperfect initial settings to improve tracking and constraint adherence, with states and inputs remaining within limits (Esfahani et al., 14 Jul 2025).

The practical guidance given for implementation is specific. Kernel selection begins with squared-exponential kernels with automatic relevance determination, hyperparameters are optimized by marginal likelihood, and Matérn kernels or additive kernels are suggested for high-dimensional $x_{k+1} = f(x_k, u_k, d_k),$73. Noise levels should be calibrated from repeated rollouts, with small jitter added for numerical stability. The EHVI reference point should be chosen slightly worse than the worst observed feasible objective vector and updated as the Pareto set improves. Initial safe parameters should emphasize conservative MPC weights, robust constraint tightenings, and large slack penalties. The gradient estimator should use batch LSTD with ridge regularization of $x_{k+1} = f(x_k, u_k, d_k),$74 and $x_{k+1} = f(x_k, u_k, d_k),$75 under near-collinearity, and KKT sensitivities should be computed at converged solutions with consistent active sets. For model mismatch, the paper recommends maintaining high slack penalties, keeping the stability objective $x_{k+1} = f(x_k, u_k, d_k),$76, and considering tube MPC parameterization in $x_{k+1} = f(x_k, u_k, d_k),$77. Objective trade-offs are described explicitly: $x_{k+1} = f(x_k, u_k, d_k),$78 versus $x_{k+1} = f(x_k, u_k, d_k),$79 trades performance against convergence to stationarity, $x_{k+1} = f(x_k, u_k, d_k),$80 stabilizes critic updates, and $x_{k+1} = f(x_k, u_k, d_k),$81 promotes closed-loop stability.

The stated limitations are also specific. Scalability is constrained by high-dimensional $x_{k+1} = f(x_k, u_k, d_k),$82, long horizons, and the cubic scaling of GP updates. MPC sensitivities are only locally valid and may become nonsmooth at active-set changes. GP-based safe exploration can be conservative. BO assumes stationary objectives, which is problematic under plant drift. Explicitly constrained multi-objective EHVI and integration with chance constraints are identified as promising directions. Another stated future direction is integration with robust identification, including learning $x_{k+1} = f(x_k, u_k, d_k),$83 or uncertainty sets in parallel with cost-based tuning. A plausible implication is that the framework is intended as a modular bridge among robust MPC design, gradient-based RL, and data-efficient black-box global search, rather than as a closed final architecture.