A Successive Constraint Approach to Solving Parameter-Dependent Linear Matrix Inequalities

TL;DR

The paper introduces a successive constraint method that efficiently solves large-scale parameter-dependent linear matrix inequalities by leveraging offline/online decomposition, enabling rapid solutions across many parameter values.

Contribution

It presents a novel approach combining offline/online decomposition to efficiently solve parameter-dependent LMIs and extends it to semidefinite programming approximations.

Findings

Significant reduction in online computation time.

Effective handling of large-scale problems.

Extension to semidefinite programming.

Abstract

We present a successive constraint approach that makes it possible to cheaply solve large-scale linear matrix inequalities for a large number of parameter values. The efficiency of our method is made possible by an offline/online decomposition of the workload. Expensive computations are performed beforehand, in the offline stage, so that the problem can be solved very cheaply in the online stage. We also extend the method to approximate solutions to semidefinite programming problems.

Equations

$$F(x;\mu ):=\sum _{q=1}^{Q_F}\left[{\theta }_q^0(\mu )+{\theta }_q^L(\mu )x\right]F_q⪰0.$$

$$\underset{x\in {\mathbb{R}}^n}{\mathrm{minimize}}c(\mu {)}^{T}x\text{subject to}F(x;\mu )⪰0.$$

$$\alpha (x;\mu ):=\underset{v\in {\mathbb{R}}^{\mathcal{N}}}{\inf }\frac{v^{T}F(x;\mu )v}{v^T{F}_{S}v},$$

$$\alpha(x;\mu)=\inf_{y\in\mathcal{Y}}\left[\theta^{0}(\mu)+\theta^{L}(\mu)x\right]^{T}y,\quad\text{where}\quad\mathcal{Y}:=\left\{y\in\mathbb{R}^{Q_{F}}\Big{|}y_{q}=\frac{v^{T}F_{q}v}{v^{T}F_{S}v},v\in\mathbb{R}^{\mathcal{N}}\right\}.$$

$${\mathcal{B}}_Q:=\prod _{i=1}^{Q_F}\left[\underset{y\in \mathcal{Y}}{\inf }{y}_q,\underset{y\in \mathcal{Y}}{\sup }{y}_q\right]\subset {\mathbb{R}}^{Q_F}.$$

$$\mathcal{Y}_{out}:=\left\{y\in\mathcal{B}_{Q}\Big{|}\left[\theta^{0}(\bar{\mu})+\theta^{L}(\bar{\mu})\bar{x}\right]^{T}y\geq\bar{\alpha}(\bar{x};\bar{\mu}),\,\forall\ (\bar{x},\bar{\mu})\in\mathcal{R}\right\}.$$

$${\alpha }_{out}(x;\mu ):=\underset{y\in {\mathcal{Y}}_{out}}{\inf }{\left[{\theta }^0(\mu )+{\theta }^L(\mu )x\right]}^{T}y,$$

$$\underset{y\in {\mathbb{R}}^{Q_F}}{\mathrm{minimize}}$$

$$A_{out}y\ge b_{out}$$

$$\underset{p\in {\mathbb{R}}^{\ell }}{\mathrm{maximize}}$$

$$A_{out}^{T}p={\theta }^0(\mu )+{\theta }^L(\mu )x$$

$$p\ge 0$$

$$\underset{p\in {\mathbb{R}}^{\ell },x\in {\mathbb{R}}^n}{\mathrm{maximize}}{b}_{out}^{T}p\text{subject to}{A}_{out}^{T}p-{\theta }^L(\mu )x={\theta }^0(\mu )\text{and}p\ge 0.$$

$${\alpha }_{in}(x;\mu ):=\underset{\bar{y}\in {\mathcal{Y}}_{in}}{\inf }{\left[{\theta }^0(\mu )+{\theta }^L(\mu )x\right]}^T\bar{y},$$

$$\underset{x\in {\mathbb{R}}^n,p\in {\mathbb{R}}^{\ell }}{\mathrm{minimize}}c(\mu {)}^{T}x\text{subject to}\left[\begin{array}{cc}{A}_{out}^T& -{\theta }^L(\mu )\end{array}\right]\left[\begin{array}{c}p\\ x\end{array}\right]={\theta }^0(\mu )\text{and}\left[\begin{array}{c}I\\ b_{out}^T\end{array}\right]p\ge \left[\begin{array}{c}0\\ 0\end{array}\right].$$

$$\underset{x\in {\mathbb{R}}^n}{\mathrm{minimize}}c(\mu {)}^{T}x\text{subject to}{\alpha }_{in}(x;\mu )\ge 0.$$

$$\underset{x\in {\mathbb{R}}^n}{\mathrm{minimize}}c(\mu {)}^{T}x\text{subject to}{\alpha }_{in}(x;\mu )>{\alpha }_{min},$$

$$\dot{y}=\Delta y+\mu y{1}_{{\Omega }_1}+1_{{\Omega }_2}u,\forall z\in \Omega ;\frac{\partial y}{\partial \eta }=0,\forall z\in \partial \Omega ;u=-x\zeta ;\zeta ={\int }_{{\Omega }_2}y.$$

$$M\dot{y}+A(x;\mu )y=0,\text{where}A(x;\mu )=A_0-\mu A_1+x{A}_2,\text{with}M,A_0,A_1,A_2\in {\mathbb{R}}^{2601\times 2601}.$$

$$F(x;\mu )⪰0⟺A(x;\mu )-\rho F_S⪰0,$$

Full text

A Successive Constraint Approach to Solving Parameter-Dependent Linear Matrix Inequalities

Robert O’Connor

Abstract:

We present a successive constraint approach that makes it possible to cheaply solve large-scale linear matrix inequalities for a large number of parameter values. The efficiency of our method is made possible by an offline/online decomposition of the workload. Expensive computations are performed beforehand, in the offline stage, so that the problem can be solved very cheaply in the online stage. We also extend the method to approximate solutions to semidefinite programming problems.

1 Introduction

Linear matrix inequalities (LMIs) are a general type of convex constraint that include linear as well as quadratic constraints [11, 12] and lead to very natural formulations of a large number of problems in control and systems theory [2, 4]. They can be solved using a wide range of existing methods [1, 5, 6, 10], but that can be expensive for large-scale problems. In particular, problems resulting from the discretization of partial differential equations can be extremely expensive due to their high dimensionality. The computations are even more expensive if parameter-dependent problems are considered and solutions are needed for a large number of parameter values. In that case traditional solution methods are extremely inefficient. We propose the construction of reduced-order models that take advantage of the parametric nature of the problem and allow us to very cheaply produce solutions for a large number of parameter values.

Let us introduce a finite-dimensional, bounded parameter domain $\mathcal{D}\in\mathbb{R}^{p}$ and the parameter-dependent LMI

[TABLE]

Here $F(x;\mu)\in\mathbb{R}^{\mathcal{N}\times\mathcal{N}}$ is a symmetric matrix that depends on both the parameter $\mu\in\mathcal{D}$ and the decision variable $x\in\mathbb{R}^{n}$ and is composed of $Q_{F}$ parameter-independent matrices $F_{q}\in\mathbb{R}^{\mathcal{N}\times\mathcal{N}}$. The parameter dependencies of $F(x;\mu)$ are given by the functions $\theta^{0}(\cdot):\mathcal{D}\rightarrow\mathbb{R}^{Q_{F}}$ and $\theta^{L}(\cdot):\mathcal{D}\rightarrow\mathbb{R}^{Q_{F}\times n}$. We use subscripts to indicate components of a vector, such that $\theta_{q}^{0}(\mu)$ is the $q$th element of $\theta^{0}(\mu)$. Similarly, we write $\theta_{q}^{L}(\mu)$ to indicate the $q$th row of $\theta^{L}(\mu)$. The symbol $\succeq$ will be used in the sense that $P\succeq 0$ indicates that the symmetric matrix $P$ is positive semi-definite.

The goal of this paper is to efficiently solve the following problems for a large number of parameter values $\mu\in\mathcal{D}$:

The strict feasibility problem: Find an $x\in\mathbb{R}^{n}$ such that $F(x;\mu)\succ 0$. 2. 2.

The semidefinite program (SDP):

[TABLE]

Rather than directly solving these problems for each new parameter value, we will solve them for only a small set of intelligently chosen parameter values. We will then use the resulting solutions to build a reduced-order model that can approximate the solution anywhere in $\mathcal{D}$.

The method that we propose can be viewed as a generalization of the successive constraint method (SCM) [3, 7], which is often used in the field of reduced basis methods to evaluate stability constants [9]. This method will allow us to very cheaply determine feasible solutions for any $\mu\in\mathcal{D}$. That is made possible by decomposing the computational workload into offline and online stages. All expensive computations will be performed in advance, during the offline stage. During the online stage the cost to solve the problem for a new parameter value will be independent of the size of the original constraint, $\mathcal{N}$. In that way the computational cost of each new solution will remain cheap even if the original constraint has very large dimensions.

Our methods are applicable to a wide range of LMIs and can also be used to extend the applicability of SCM. In the context of reduced basis methods, applications could involve bounding stability constants with respect to parameter-dependent norms or the selection of Lyapunov functions for the computation of error bounds [8]. In Section 4 we present an example in which we optimize a system while ensuring that it remains stable.

2 Reduced-order modeling for strict feasibility

SCM was originally designed to approximate coercivity constants. We will apply a modified version of SCM to the coercivity constant

[TABLE]

where $F_{S}\in\mathbb{R}^{\mathcal{N}\times\mathcal{N}}$ is a fixed symmetric positive-definite matrix. From the definition it is clear that $\alpha(x;\mu)\geq 0$ is equivalent to $F(x;\mu)\succeq 0$ for all symmetric positive definite matrices $F_{S}$. Nevertheless, an appropriate choice of $F_{S}$ could be beneficial from a numerical point of view. If we are dealing with PDE discretizations, it can be advantageous to choose a matrix associated with an energy norm.

The first step in applying SCM is reformulating the coercivity constant as follows:

[TABLE]

This formulation has the advantage that the complexity of the problem has been shifted to the definition of the set $\mathcal{Y}$. That allows us to compute lower and upper bounds for $\alpha(x;\mu)$ by approximating $\mathcal{Y}$.

A lower bound for $\alpha(x;\mu)$ can be derived by approximating $\mathcal{Y}$ from the outside. A bounded but primitive approximation for $\mathcal{Y}$ is given by

[TABLE]

$\mathcal{B}_{Q}$ will generally be much larger than $\mathcal{Y}$ so we will add constraints to restrict it and better approximate $\mathcal{Y}$. Let us assume that we have computed lower bounds $\bar{\alpha}(\bar{x};\bar{\mu})\leq\alpha(\bar{x};\bar{\mu})$ for each pair $(\bar{x},\bar{\mu})$ in a predetermined set $\mathcal{R}\subset\mathbb{R}^{n}\times\mathcal{D}$. The set $\mathcal{Y}$ is then contained in the set

[TABLE]

Using $\mathcal{Y}_{out}$ we can define

[TABLE]

such that $\alpha_{out}(x;\mu)\leq\alpha(x;\mu)$ for any $x\in\mathbb{R}^{n}$ and any $\mu\in\mathcal{D}$.

The lower bound $\alpha_{out}(x;\mu)$ can be formulated as the solution to a linear programming problem. To see that we note that $\mathcal{Y}_{out}$ is a linearly constrained subset of $\mathbb{R}^{Q_{F}}$. As a result we can find a matrix $A_{out}\in\mathbb{R}^{\ell\times Q_{F}}$ and a vector $b_{out}\in\mathbb{R}^{\ell}$ such that $\mathcal{Y}_{out}=\{y\in\mathbb{R}^{Q_{F}}|A_{out}y\geq b_{out}\}$. Here the symbol $\geq$ indicates a componentwise comparison of vectors and $\ell\ll\mathcal{N}$ is the number of constraints. The lower bound $\alpha_{out}(x;\mu)$ can be calculated, for any $x\in\mathbb{R}^{n}$ and any $\mu\in\mathcal{D}$, as the solution to the linear program (P) or its dual (D):

$\displaystyle\underset{y\in\mathbb{R}^{Q_{F}}}{\text{minimize}}$

$\displaystyle\left[\theta^{0}(\mu)+\theta^{L}(\mu)x\right]^{T}y$

(P)

subject to

$\displaystyle A_{out}y\geq b_{out}$

$\displaystyle\underset{p\in\mathbb{R}^{\ell}}{\text{maximize}}$

$\displaystyle b_{out}^{T}p$

(D)

subject to

$\displaystyle A_{out}^{T}p=\theta^{0}(\mu)+\theta^{L}(\mu)x$

$\displaystyle p\geq 0$

To find an optimal value of $x$ we will maximize (D) over all $x\in\mathbb{R}^{n}$ using the problem

[TABLE]

If $b_{out}^{T}p$ is strictly greater than [math], then the associated $x$ is guaranteed to strictly satisfy (L). Even if that is not the case, (RF) is always feasible. To show that we note that the boundedness of $\mathcal{Y}_{out}\subset\mathcal{B}_{Q}$ implies the boundedness of (P). By duality (D) is then feasible for any $x\in\mathbb{R}^{n}$ and any $\mu\in\mathcal{D}$ which implies that (RF) is feasible for any $\mu\in\mathcal{D}$. (RF) is also bounded if $\max_{x\in\mathbb{R}^{n}}\alpha(x;\mu)$ is bounded; otherwise, we will be content with any $x$ associated with a large, positive value of $b_{out}^{T}p$. The main reason to use (RF) is that it is cheap to solve: It is a small linear program with $n+\ell$ variables, $Q_{F}$ equality constraints, and $\ell$ inequality constraints.

Since it is not possible to work directly with the infinite set $\mathcal{D}$, we introduce a large set $\Xi\subset\mathcal{D}$ of discrete points that are representative of $\mathcal{D}$ and a much smaller set $\mathcal{C}_{k}\subset\Xi$ with cardinality $k$. The set $\mathcal{R}$ will be made up of parameter values $\bar{\mu}\in\Xi$ and associated approximate solutions $\bar{x}\in\mathbb{R}^{n}$. For parameter values $\bar{\mu}$ that are in $\mathcal{C}_{k}$ we will use solutions $\bar{x}$ to the original problem (L) and set $\bar{\alpha}(\bar{x};\bar{\mu})=\alpha(\bar{x};\bar{\mu})$. For parameter values in the much larger set $\Xi\setminus\mathcal{C}_{k}$ we will use solutions $\bar{x}$ to (RF), and set $\bar{\alpha}(\bar{x};\bar{\mu})=\alpha_{out}(\bar{x};\bar{\mu})$. The pairs $(\bar{x},\bar{\mu})$ that will be included in $\mathcal{R}$ will depend on the parameter $\mu$ for which we would like a feasible solution. For two natural numbers $M_{C}$ and $M_{\Xi}$ we define $\mathcal{R}$ to be a set of $M_{C}$ pairs $(\bar{x},\bar{\mu})$ with $\bar{\mu}\in\mathcal{C}_{k}$ and $M_{\Xi}$ pairs $(\bar{x},\bar{\mu})$ with $\bar{\mu}\in\Xi\setminus\mathcal{C}_{k}$. In both cases we will choose the pairs $(\bar{x},\bar{\mu})$ with $\bar{\mu}$ closest to $\mu$ in a predetermined metric. If $M_{C}>k$, we simply include all of $\mathcal{C}_{k}$ in $\mathcal{R}$.

To build our model we will make use of the greedy algorithm that was introduced in the context of reduced basis methods [13] and plays a vital role in SCM. The algorithm is initiated by choosing a small initial set $\mathcal{C}_{k}$. For each $\bar{\mu}\in\mathcal{C}_{k}$ the problem (L) is solved once to determine the associated values of $\bar{x}$ and $\bar{\alpha}(\bar{x};\bar{\mu})$. For $\bar{\mu}\in\Xi\setminus\mathcal{C}_{k}$ we initially set $\bar{x}=0$ and $\bar{\alpha}(\bar{x};\bar{\mu})=-\infty$. The model is then improved iteratively. In each iteration we will solve (RF) for each $\bar{\mu}\in\Xi\setminus\mathcal{C}_{k}$ and update the stored values of $\bar{x}$ and $\bar{\alpha}(\bar{x};\bar{\mu})$. The $\bar{\mu}$ associated with the smallest value of $\bar{\alpha}(\bar{x};\bar{\mu})$ is then added to $\mathcal{C}_{k}$, and (L) is solved to update the stored values of $\bar{x}$ and $\bar{\alpha}(\bar{x};\bar{\mu})$. The process terminates when the smallest value of $\bar{\alpha}(\bar{x};\bar{\mu})$ is larger than some positive tolerance.

Solving problems with our method involves two stages. During the offline stage the greedy algorithm builds the model. The expensive operations in the offline stage are $k$ solutions of (L) and nearly $k|\Xi|$ solutions of (RF). Here $k$ is the cardinality of the final $\mathcal{C}_{k}$ and $|\Xi|$ is the cardinality of $\Xi$. During the online stage (RF) can be constructed and solved very cheaply for any new parameter value. The online computational cost is independent of $\mathcal{N}$ to ensure that it remains cheap even if $\mathcal{N}$ is very large.

During the greedy algorithm it is necessary to solve (L) for various parameter values. That can be done using a variety of methods including algorithms that solve semidefinite programming problems [1, 5, 6, 10]. In this section we propose a method that is based on SCM and allows us to reuse information that we have already computed.

For a predetermined, finite subset $\mathcal{Y}_{in}$ of $\mathcal{Y}$ we define

[TABLE]

such that $\alpha_{in}(x;\mu)\geq\alpha(x;\mu)$ for all $x\in\mathbb{R}^{n}$ and $\mu\in\mathcal{D}$. Here we will make use of a variation of the greedy algorithm. In each iteration we solve $\max_{x\in\mathbb{R}^{n}}\alpha_{in}(x;\mu)$ to find the optimal $x$. For that value of $x$ we then compute $\alpha(x;\mu)$ and an element $y$ of $\mathcal{Y}$ for which the infimum in (2) is reached. That $y$ is then added to $\mathcal{Y}_{in}$ to improve the upper bound $\alpha_{in}(\cdot;\mu)$. The algorithm terminates when $\alpha(x;\mu)$ is sufficiently large or sufficiently close to $\alpha_{in}(x;\mu)$.

In each iteration of the algorithm both the small linear program associated with $\max_{x\in\mathbb{R}^{n}}\alpha_{in}(x;\mu)$ and the $\mathcal{N}$-dimensional eigenvalue problem associated with $\alpha(x;\mu)$ need to be solved once. The dominant cost is that of the eigenvalue problems. To reduce the number of iterations and eigenvalue solves we will store the last value of $\bar{y}$ that is calculated each time we solve (L) for a new $\mu\in\mathcal{C}_{k}$. The set $\mathcal{Y}_{in}$ can then be initialized using those values of $\bar{y}$.

3 The semidefinite programming problem

In this section we will show how our methods can be used to approximate solutions to the semidefinite program (S). The idea is to minimize $c(\mu)^{T}x$ over all $x\in\mathbb{R}^{n}$ that satisfy $\alpha_{out}(x;\mu)\geq 0$. That is done using the following problem:

[TABLE]

Here $A_{out}$, $b_{out}^{T}$ and $\alpha_{out}(x;\mu)$ are all the same as in Section 2 except that the stored values of $\bar{x}$ will be calculated using (S) and (RS) rather than (L) and (RF). The new constraint $b_{out}^{T}p\geq 0$ guarantees that $\alpha_{out}(x;\mu)\geq 0$.

Let us write $\mathcal{J}_{out}$ and $\mathcal{J}$ to denote the optimal values of (RS) and (S), respectively. Based on the fact that the optimal $x\in\mathbb{R}^{n}$ for (RS) will be feasible for (S) we know that $\mathcal{J}_{out}\geq\mathcal{J}$. Ideally $\mathcal{J}_{out}$ will also be a good approximation to $\mathcal{J}$ despite being much cheaper to compute: It requires the solution of a linear program with $n+\ell$ variables, $Q_{F}$ equality constraints, and $\ell+1$ inequality constraints.

For any $\mu\in\mathcal{D}$, the boundedness of (S) implies the boundedness of (RS), but the question of feasibility is more complicated. To ensure feasibility we first build a model (RF) to find strictly feasible solutions. We then convert that model to the form of (RS). In doing so we consider the stored values of $\bar{x}\in\mathbb{R}^{n}$, which should be feasible, to be approximate solutions to (S) and we reinitialize $\mathcal{C}_{k}$ as an empty set. We can then improve the model using another greedy algorithm.

The goal of this greedy algorithm is to reduce the error $\mathcal{J}_{out}-\mathcal{J}$. Since it is too expensive to compute $\mathcal{J}$ online, we will compute a lower bound for it by considering the problem

[TABLE]

Here $\alpha_{in}(x;\mu)$ is defined like in (6) with $\mathcal{Y}_{in}$ initially being the set of $\bar{y}\in\mathcal{Y}$ associated with parameter values in $\mathcal{C}_{k}$. The optimal value of (ER), which we will denote $\mathcal{J}_{in}$, is by definition a lower bound for $\mathcal{J}$. That allows us to bound the error $\mathcal{J}_{out}-\mathcal{J}\geq 0$ from above using $\mathcal{J}_{out}-\mathcal{J}_{in}$. Like the previous approximations, this bound can be computed online with a number of computations that is independent of $\mathcal{N}$. Although (ER) is always feasible, it could be unbounded. That can be fixed by increasing $\mathcal{C}_{k}$ and consequently $\mathcal{Y}_{in}$.

During each iteration of the greedy algorithm we will solve both (RS) and (ER), compute the difference $\mathcal{J}_{out}-\mathcal{J}_{in}$, and update both the estimates $\bar{x}\in\mathbb{R}^{n}$ and the associated values $\bar{\alpha}(\bar{x};\bar{\mu})$ for each $\bar{\mu}\in\Xi\setminus\mathcal{C}_{k}$. The $\bar{\mu}$ that produced the maximum value of $\mathcal{J}_{out}-\mathcal{J}_{in}$ is then added to $\mathcal{C}_{k}$ and a more accurate solution to (S) is computed. That allows us to update both (RS) and (ER). This process terminates when the maximum value of $\mathcal{J}_{out}-\mathcal{J}_{in}$ is below a desired tolerance.

The major computational burdens in the offline stage of our method are the construction of the model (RF), $k$ solutions of (S), and nearly $k|\Xi|$ solutions of both (RS) and (ER). During the online stage an approximate solution can be determined for any given $\mu\in\mathcal{D}$ by simply constructing and solving (RS). The error can also be bounded online by solving (ER).

The more accurate solutions that we need can be computed using either general solvers for semidefinite programming problems [1, 5, 6, 10] or a method that is similar to the one presented in Section 2. For the latter we will consider the following problem with a fixed value of $\mu\in\mathcal{D}$:

[TABLE]

where $\alpha_{min}>0$ is some small constant that is used as a tolerance. We again use a sort of greedy algorithm. In each iteration we solve (T) to get the approximate solution $x$ and update $\mathcal{Y}_{in}$ by adding the new value of $y$. The process is repeated until the stopping condition $\alpha(x;\mu)>0$ is satisfied. If $\alpha_{min}$ is sufficiently small, the resulting $x$ should be a good approximation to the true solution.

4 Numerical example: System stabilization

To test our method we will consider a reaction-diffusion equation on the unit square, $\Omega$, depicted in Figure 1(a). To facilitate the explanation we consider a problem with just one parameter $\mu\in\mathcal{D}:=[0,3]$ and one decision variable $x\in\mathbb{R}$. The problem is related to the following system:

[TABLE]

Here $\mathbf{1}_{S}$ indicates the characteristic function of a set $S$ and $\eta$ is the outward-pointing unit vector on the boundary $\partial\Omega$. The decision variable $x$ determines the feedback gain from the system’s output $\zeta(t)$ to its input $u(t)$ as depicted in Figure 1(b). We use linear finite elements to discretize the problem in space. That gives us the following $2601$-dimensional closed-loop semi-discrete system

[TABLE]

Here $M$ is symmetric positive definite and $A(x;\mu)$ is symmetric. Under these assumptions the system (8) is strictly stable for a given $\mu\in\mathcal{D}$ and $x\in\mathbb{R}$ iff $A(x;\mu)\succ 0$. Our goal will be to minimize the cost of the control, by keeping $x$ small, while ensuring strict stability. Noting that $A_{0}$, which is associated with the operator $-\Delta$, is positive semidefinite and that $A_{0}+A_{1}$ is positive definite, we define $F(x;\mu):=(1-\rho)A_{0}+(-\mu-\rho)A_{1}+xA_{2}$ and $F_{S}:=A_{0}+A_{1}\succ 0$ for some small $\rho>0$. It then holds that

[TABLE]

and hence the system in (8) is strictly stable if $F(x;\mu)\succeq 0$. We will search for a minimal stabilizing gain $x\in\mathbb{R}$ as the solution to (S).

For our model we will use values of $M_{C}=4$ and $M_{\Xi}=3$. $\mathcal{C}_{k}$ is initialized using the smallest and the largest values in $\mathcal{D}$. For this particular problem that guarantees that (RS) is feasible for all $\mu\in\mathcal{D}$. As a result we can start directly with (RS) and do not need to build the model given in (RF). Another modification of our method is that we will work with the relative error given by $(\mathcal{J}_{out}-\mathcal{J}_{in})/\mathcal{J}_{in}$ rather than the error $\mathcal{J}_{out}-\mathcal{J}_{in}$. That is possible because $\mathcal{J}_{in}$ will always be strictly positive.

Figure 2(a) shows the convergence of the solver described in Section 3 for the computation of more accurate solutions to (S). Here the error is measured as $\alpha_{in}(\bar{x};\bar{\mu})-\alpha(\bar{x};\bar{\mu}).$ Two examples are shown as well as the worst case over $30$ random parameter values. We recall that the most expensive part of this method is the eigenvalue solves and that one such solve is needed for each iteration. Figure 2(b) shows the convergence of the greedy algorithm for two different parameter domains $\mathcal{D}$. The plotted values are the worst relative errors over the respective sets $\Xi$. For $\mathcal{D}=[0,3]$ we set $\Xi$ to be $300$ uniformly distributed points and for $\mathcal{D}=[0,1.5]$ we use $150$ points. Figure 2(c) compares the results when $\mathcal{C}_{k}$ is constructed using the greedy algorithm or simply as uniformly distributed points. The plot shows the worst relative errors over $100$ random values of $\mu\in\mathcal{D}$. The greedy algorithm performs slightly better than uniform distributions when $\mathcal{C}_{k}$ is large and also has the advantage that it is iterative.