Zeroth-order optimisation on subsets of symmetric matrices with application to MPC tuning

Let $X^{*}\in\mathbb{Q}^{n}$ be a stationary point of (1), and $f^{*}\triangleq f(X^{*})$. In addition, define $U_{0:k}\triangleq\mbox{diag}\{U_{0},\dots,U_{k}\}$, a random matrix composed by i.i.d. variables $\{U_{k}\}_{k\in\mathbb{N}_{0}}$ associated with each iteration of the scheme.

We are now in a position to state the main result.

Proof: See Appendix B. $\blacksquare$

A corollary from Theorem 1 can be obtained which provides expressions for $\mu$, $h_{k}$, and $N$ that ensure a given level of accuracy for the ZO-RMS algorithm.

Proof: Note that

Suppose the number of steps $N$ is fixed, we can choose $\mu$ and $h_{k}$ as in (5) and further obtain the following bound

Therefore, in order to satisfy $\mathbb{E}_{U_{0:N-1}}\{f(\hat{X}_{N})\}-f^{*}\leq\epsilon$, we need $N\geq\frac{L_{0}^{2}(f)\bar{r}^{2}}{\epsilon^{2}}(n^{4}+2n^{3}+5n^{2}+4n)-1$, concluding the proof. $\blacksquare$

The complexity bound (6) in Corollary 1 corresponds to the number of iterations in which the ZO-RMS algorithm guarantees a given accuracy of $\epsilon$, provided the step size $h_{k}$ and parameter $\mu$ are chosen as per (5). Certainly, the expressions for $h_{k},\mu$, and $N$ depend on the Lipschitz constant $L_{0}(f)$ and $\bar{r}$ which may be hard to obtain explicitly depending on the application. However, these can be numerically bounded as we explain further below in Section 6.1. Since the oracle (3) is random, the guarantees of the ZO-RMS algorithm hold in expectation. In essence, Corollary 1 provides sufficient conditions on the oracle’s precision $\mu$ and step size $h_{k}$ such that we can get sufficiently close to the optimal cost function $f^{*}$ in $N$ iterations, in average.

We now compare our complexity bound (6) to the bounds in [27] which are applicable to (1) (for convex $f$) after grouping the elements of the tuning matrix $X$ into a single vector. Particularly, since $X$ is symmetric, we group the lower triangular elements into a vector of size $n(n+1)/2$ (or often called half vectorization), therefore the results in [27] apply and the following complexity bound holds (cf. Theorem 6 in [27])

Fig. 1 depicts (6) and (7) for different values of $n$. We can see that our bound (6) is less conservative than (7) for all $n\in\mathbb{N}$, and tends to (7) as $n\rightarrow\infty$. We emphasise that our optimisation framework is directly tailored to matrix parameters, and the reduction in conservatism with respect to a vector approach such as [27] follows from exploiting the special structure of the underlying matrix space, namely the symmetry, and the use of iterates that are tailored to it. The reduction in conservativeness is significant, for instance, for $3\times 3$ matrices (i.e. $n=3$), $N_{\text{\tiny\cite[cite]{[\@@bibref{}{nesspo17}{}{}]}}}/N_{\text{\tiny Cor.\ref{coro:complexity-bound}}}=2.08$. This translates in a sufficient condition that guarantees a level of accuracy in twice less iterations.

Compared to the convex optimisation problem considered in Section 3.1, the analysis for constrained non-convex problems is more challenging. Particularly, constrained convex problems typically focus on the optimality gap $\mathbb{E}\left\{f(x)\right\}-f^{*}$ to measure the convergence rate (as we do in Corollary 1). On the other hand, Nesterov in [27] provided complexity bounds for unconstrained non-convex problems using zeroth order iterates, where the optimisation variable is a vector in a subset of $\mathbb{R}^{n}$. In these unconstrained non-convex problems, the object $\mathbb{E}\left\{\left\|\nabla f(x)\right\|^{2}\right\}$ is the typical measure for stationarity. We emphasise that the bounds for unconstrained non-convex problems in [27] are not applicable to our constrained problem (1), since we need to search over the matrix space $\mathbb{Q}^{n}$ and use iterates of the form (2) that project onto $\mathbb{Q}^{n}$. For constrained non-convex problems like ours, a fitting alternative to $\mathbb{E}\left\{\left\|\nabla f(x)\right\|^{2}\right\}$ is to consider the so-called gradient mapping, see e.g. [21, 13], which is defined as follows

where $\mathcal{O}_{\mu}$ is the random zeroth order oracle in (3), and recall that $\pi_{\mathbb{Q}^{n}}$ denotes the Euclidean projection onto the closed convex set $\mathbb{Q}^{n}$. The natural interpretation of $P_{\mathbb{Q}^{n}}$ is that it represents the projected gradient, which offers a feasible update from the previous point $X$. The main goal in this section is to provide complexity bounds for the ZO-RMS algorithm (2) in terms of bounding $\mathbb{E}\left\{\left\|P_{\mathbb{Q}^{n}}(X,\mathcal{O}_{\mu},h)\right\|_{F}^{2}\right\}$ when solving the constrained optimisation problem (1) with $f$ non-convex.

Before stating our results, we impose an additional assumption.

Assumption 2 essentially bounds the variance of the random oracle $\mathcal{O}_{\mu}$. This assumption is often adopted in non-convex problems, see e.g. [13]. If constructing $\sigma$ explicitly is of interest, we can construct it as follows. We note that $\mathcal{O}_{\mu}$ is an unbiased estimator of $\nabla f_{\mu}(X)$ (cf. (18)), i.e. $\mathbb{E}\left\{\mathcal{O}_{\mu}(X,U)\right\}=\nabla f_{\mu}(X)$, then $\mathbb{E}\left\{\left\|\mathcal{O}_{\mu}(X,U)-\nabla f_{\mu}(X)\right\|_{F}^{2}\right\}\leq\mathbb{E}\left\{\left\|\mathcal{O}_{\mu}(X,U)\right\|_{F}^{2}\right\}\stackrel{{\scriptstyle\text{\tiny Thm.\ref{theo:4-matrix}}}}{{\leq}}(1/4)L_{0}^{2}(f)(n^{4}+2n^{3}+5n^{2}+4n)\triangleq\sigma^{2}$.

The main result for the non-convex case is stated below. It can be seen as the non-convex counterpart of Theorem 1.

Proof: See Appendix C. $\blacksquare$

The next corollary provides a choice of $\mu,h_{k}$, and $N$ that ensure a given level of accuracy for the ZO-RMS algorithm up to an error of order $\sigma^{2}$.

Proof: Recall that $f_{\mu}$ is a smooth approximation of $f$, and the gap in this approximation can be upper bounded by $|f_{\mu}(X)-f(X)|\leq\mu L_{0}(f)\sqrt{\frac{n^{2}+n}{2}}$, see Lemma 2-(III) in the appendix. Then, to bound this gap by $\varepsilon$ we need to choose $\mu$ as per (9). For this choice of $\mu$, which we denote $\bar{\mu}$, $C(\bar{\mu})=\frac{L_{0}^{4}(f)}{4\varepsilon}\left(\frac{n^{2}+n}{2}\right)(n^{4}+2n^{3}+5n^{2}+4n)$. For a constant step size, the right-hand size of (8) becomes

Let $\rho(h)\triangleq\frac{L_{0}(f)\bar{r}}{(N+1)h}+C(\bar{\mu})h$. Minimising $\rho(h)$ in $h$ leads to $h^{*}=\sqrt{L_{0}(f)\bar{r}/[C(\bar{\mu})(N+1)]}$, which corresponds to (10). For this choice of step size, we get $\rho(h^{*})=\sqrt{4L_{0}(f)\bar{r}C(\bar{\mu})/(N+1)}$. Note that $\mathbb{E}\{\left\|P_{\mathbb{Q}^{n}}(X_{D},\mathcal{O}_{\mu}(X_{D},U_{D}),h)\right\|_{F}^{2}\}\leq\frac{1}{N+1}\sum_{k=0}^{N}\mathbb{E}_{U_{0:k}}\{\left\|P_{\mathbb{Q}^{n}}(X_{k},\mathcal{O}_{\mu}(X_{k},U_{k}),h)\right\|_{F}^{2}\}\leq\rho(h^{*})+\sigma^{2}$. Then, we can guarantee that $\rho(h^{*})\leq\delta$ in $4L_{0}(f)\bar{r}C(\bar{\mu})/\delta^{2}$ iterations, which corresponds to (11), completing the proof. $\blacksquare$

Our analysis shows that ZO-RMS for constrained non-convex problems can suffer an additional error of order $\sigma^{2}$ which does not appear in the convex case, see Corollary 1. This effect is consistent with the literature on constrained non-convex problems, where this effect has been reported for different algorithms, see e.g. [21, 13, 20].

A schematic representation of the diesel air-path is shown in Fig. 2. The dynamics of a diesel engine air-path are highly non-linear, see e.g. [41], and can be captured in the general form

where $x_{k}\in\mathbb{R}^{n_{x}},y_{k}\in\mathbb{R}^{n_{y}}$, and $u_{k}\in\mathbb{R}^{n_{u}}$ are the state, output, and input of the system respectively, at time instant $k\in\mathbb{N}_{0}$. The engine operational space is typically parametrised by $\theta\triangleq(\omega_{e},\dot{m}_{f})\in\Theta\subseteq\mathbb{R}^{2}$, where $\omega_{e}$ denotes the engine speed, $\dot{m}_{f}$ denotes the fuel rate, and $\Theta$ is the engine operating space defined as $\Theta\triangleq\{\theta\in\mathbb{R}^{2}:\omega_{e,\min}\leq\omega_{e}\leq\omega_{e,\max},\dot{m}_{f,\min}\leq\dot{m}_{f}\leq\dot{m}_{f,\max}\}$, for some $\omega_{e,\min},\omega_{e,\max},\dot{m}_{f,\min},\dot{m}_{f,\max}\in\mathbb{R}$. Given these highly non-linear dynamics, a common approach to control-oriented modelling of the air-path is to generate linear models trimmed at various operating points $\theta$. Particularly, we follow [30] and use twelve models to approximate the operating space uniformly. That is, the engine operating space $\Theta$ is divided into twelve regions as per Fig. 3, with a linear model representing the engine dynamics in each region. For commercial in confidence purposes, we only show normalised axes in every plot. We emphasise that these regions are chosen to provide adequate coverage over the range of operating points encountered along a drive cycle. The control-oriented model state consists in the intake manifold pressure $p_{\text{im}}$ (also known as boost pressure), the exhaust manifold pressure $p_{\text{em}}$, the compressor flow $W_{\text{comp}}$, and the EGR rate $y_{\text{EGR}}$. The input consists of the throttle position $u_{\text{thr}}$, EGR valve position $u_{\text{EGR}}$, and VGT valve position $u_{\text{VGT}}$. Lastly, the measured output is $(p_{\text{im}},y_{\text{EGR}})$.

For a given operating point $(\omega_{e}^{\sigma},\dot{m}_{f}^{\sigma})$, $\sigma\in\{1,\dots,12\}$, the engine control unit (ECU) applies certain steady-state actuator values. We denote by $\bar{u}^{\sigma}\in\mathbb{R}^{3},\bar{x}^{\sigma}\in\mathbb{R}^{4}$, and $\bar{y}^{\sigma}\in\mathbb{R}^{2}$ the steady-state values of the input, state, and output respectively at each operating condition $(\omega_{e}^{\sigma},\dot{m}_{f}^{\sigma})$. Then, by following the system identification procedure detailed in [35], a linear representation of the non-linear diesel air-path (12) trimmed around the grid point $(\omega_{e}^{\sigma},\dot{m}_{f}^{\sigma})$ is given by

where $\sigma\in\{1,\dots,12\}$, $\tilde{x}_{k}=x_{k}-\bar{x}^{\sigma}\in\mathbb{R}^{4}$ is the perturbed state around the corresponding steady-state $\bar{x}^{\sigma}$, $\tilde{u}_{k}=u_{k}-\bar{u}^{\sigma}\in\mathbb{R}^{3}$ is the perturbed input around the corresponding steady-state input $\bar{u}^{\sigma}$, and $\tilde{y}_{k}=y_{k}-\bar{y}^{\sigma}\in\mathbb{R}^{2}$ is the perturbed output around the corresponding steady-state output $\bar{y}^{\sigma}$.

For each operating point $(\omega_{e}^{\sigma},\dot{m}_{f}^{\sigma})$, $\sigma\in\{1,\dots,12\}$, and corresponding model (13) associated to that region, we formulate an MPC with augmented integrator (see e.g. [29, 42]). To this end, define the augmented state $\mathbf{x}_{k}=(\tilde{x}_{k},e_{k})$, where $e_{k}$ is the integrator state with dynamics $e_{k+1}=-C^{\sigma}\tilde{x}_{k}+e_{k}$. Define the cost function as

where $H\in\mathbb{N}$ is the prediction horizon, $\mathbf{u}\triangleq\{\tilde{u}_{k},\tilde{u}_{k+1},\dots,\tilde{u}_{k+H-1}\}$ is the sequence of control values applied over the horizon $H$, and ${P}^{\sigma}\in\mathbb{S}^{6},{Q}^{\sigma}\in\mathbb{S}^{6}$, and $R^{\sigma}\in\mathbb{S}^{3}$ are real symmetric matrices containing the tuning weights. Particularly, we further assume that $R^{\sigma}$ is positive definite, and that $P^{\sigma}$ and $Q^{\sigma}$ are positive semidefinite.

The corresponding MPC problem is stated below,

where $S_{x}\in\mathbb{R}^{18\times 6}$, $S_{u}\in\mathbb{R}^{18\times 3}$, $b_{i}\in\mathbb{R}^{18}$, $S_{H}\in\mathbb{R}^{12\times 6}$, and $b_{H}\in\mathbb{R}^{12}$ are given by the relevant state and input constraints. The solution to the above optimisation problem yields the optimising control sequence $\mathbf{u}^{*}=\{\tilde{u}^{*}_{k},\tilde{u}^{*}_{k+1},\dots,\tilde{u}^{*}_{k+H-1}\}$. The first element of the sequence $\mathbf{u^{*}}$ is applied to the system and the whole process is repeated as $k$ is incremented.

For transient operations between operating points, we use a switching LTI-MPC architecture so that the controllers are switched based on the current operating condition. As in [31], the switching LTI-MPC strategy selects the MPC controller at the nearest operating point to the current operating condition.

Under the above scenario, our goal is to use the ZO-RMS algorithm to tune the MPC weighting matrices $\{P^{\sigma},Q^{\sigma},R^{\sigma}\}$ in order to achieve a satisfactory tracking performance for a given prediction horizon $H\in\mathbb{N}$. Since there are twelve available controllers to tune, we could potentially use the ZO-RMS algorithm to tune all of them. However, since experiments are costly, we will focus on tuning only the controllers that have poor performance given an initial choice of tuning parameters. That is, we utilise an initial calibration denoted by $\{P_{0}^{\sigma},Q_{0}^{\sigma},R_{0}^{\sigma}\}$ for every controller $\sigma\in\{1,\dots,12\}$, which we call baseline controller. The choice of $\{P_{0}^{\sigma},Q_{0}^{\sigma},R_{0}^{\sigma}\}$ may be based on model dynamics, experience, or using ad-hoc guidelines as per [12]. Then, by looking at a drive cycle response with the baseline controller, we detect the controllers that have poor performance. Let us denote by $\sigma^{*}$ a controller with unsatisfactory performance and that we attempt to tune with the ZO-RMS algorithm. Therefore, in this context, we want to solve

where $g$ is the tracking performance of the MPC defined in (15) below, which depends on a closed-loop response of interest $y$, and $\mathbb{S}^{4}_{+}\subset\mathbb{S}^{4}$ and $\mathbb{S}^{3}_{++}\subset\mathbb{S}^{3}$ denote the set of symmetric positive semi-definite matrices and the set of symmetric positive definite matrices, respectively.

As mentioned above, we use the tracking error as the measure of performance for the MPC controller, that is,

where $M\in\mathbb{N}$ is the experiment length, $y_{k}^{\text{ref}}$ is the vector containing the boost pressure and EGR rate references, respectively, and $y_{k}(P^{\sigma^{*}},Q^{\sigma^{*}},R^{\sigma^{*}})$ is the process measured output when the $\sigma^{*}$-th controller is using tuning parameters $P^{\sigma^{*}},Q^{\sigma^{*}},R^{\sigma^{*}}$, and the rest of controllers are using the baseline parameters $\{P_{0},Q_{0},R_{0}\}$. We emphasise that for the simulations and experiments below we input normalised signals in (15) so that they are weighted equally.

We now show how to implement the ZO-RMS algorithm for the MPC tuning problem (14). First note that (14) fits the general optimisation problem (1) with $X=\mbox{diag}\{P^{\sigma^{*}},Q^{\sigma^{*}},R^{\sigma^{*}}\}$, $f(X)\triangleq g(y(X))$, and

Therefore, we can apply ZO-RMS to solve (14). The overall implementation of ZO-RMS in this context is depicted in Fig. 4. the ZO-RMS algorithm iteratively updates the MPC parameters $\{P^{\sigma^{*}},Q^{\sigma^{*}},R^{\sigma^{*}}\}$ to minimise $f(P^{\sigma^{*}},Q^{\sigma^{*}},R^{\sigma^{*}})$, which is calculated from closed-loop response experiments carried out between the MPC and diesel engine.

It remains to show how to compute the oracle and projection in (2). Particularly, the oracle in this context is computed as per (3) with $U=\mbox{diag}\{U^{P},U^{Q},U^{R}\}$, where $U^{P},U^{Q}\in\mathbb{G}^{4}$ and $U^{R}\in\mathbb{G}^{3}$ (see Definition 1). We emphasise that, at each iteration step, the oracle is computed once the entire closed-loop response is available. Particularly, two experiments are required to compute the oracle, one with parameters $\{P_{k}^{\sigma^{*}}+\mu U_{k}^{P},Q_{k}^{\sigma^{*}}+\mu U_{k}^{Q},R_{k}^{\sigma^{*}}+\mu U_{k}^{R}\}$ and one with $\{P_{k}^{\sigma^{*}},Q_{k}^{\sigma^{*}},R_{k}^{\sigma^{*}}\}$. Once the two closed-loop responses are finished, then ZO-RMS computes the oracle and next update $\{P_{k+1}^{\sigma^{*}},Q_{k+1}^{\sigma^{*}},R_{k+1}^{\sigma^{*}}\}$. Two new closed-loop experiments are then performed with these new controller parameters and the process continues iteratively in this fashion.

Since $X$ is block diagonal, the projection onto $\mathbb{Q}^{n}$ is computed as $\pi_{\mathbb{Q}^{n}}\{X\}=\mbox{diag}\{\pi_{\mathbb{S}^{4}_{+}}\{P\},\pi_{\mathbb{S}^{4}_{+}}\{Q\},\pi_{\mathbb{S}^{3}_{++}}\{R\}\}$, where $\pi_{\mathbb{S}^{4}_{+}}$ and $\pi_{\mathbb{S}^{3}_{++}}$ denote the Euclidean projection onto $\mathbb{S}^{4}_{+}$ and $\mathbb{S}^{3}_{++}$, respectively. To compute them, we follow the well known results in [15]. That is, let $K=V\Lambda V^{\top}$ be the eigenvalue decomposition of a matrix $K$. Then, $\pi_{\mathbb{S}^{+}}\{K\}\triangleq V\max\{0,\Lambda\}V^{\top}$ and $\pi_{\mathbb{S}^{++}}\{K\}\triangleq V\max\{d,\Lambda\}V^{\top}$, $d>0$.