Conic Optimization

To solve polynomial sum-of-squares optimization problems via semidefinite programming, CaΣoS provides both a high-level and low-level interface for convex optimization. These interfaces extend CasADi's qpsol and conic syntax. Several conic solvers are supported (but, unlike CasADi, the solvers must be installed separately and be accessible through MATLAB).

Note

SeDuMi does not support quadratic cost functions.

High-level interface

The high-level interface solves convex problems of the form

$$\begin{array}{l c r} \min & f(x,p), & x = (x_\mathrm{l}, x_\mathrm{c}) \\\ \text{s.t.} & g_\mathrm{lb} \leq g_\mathrm{l}(x,p) \leq g_\mathrm{ub}, & g_\mathrm{c}(x,p) \succeq_{\mathcal K_c} g_\mathrm{cb} \\\ \text{and} & x_\mathrm{lb} \leq x_\mathrm{l} \leq x_\mathrm{ub}, & x_\mathrm{c} \succeq_{\mathcal{K}_x} x_\mathrm{cb} \end{array}$$

where $f$ is a convex quadratic function in $x$, the constraints $g_\mathrm{l}$ and $g_\mathrm{c}$ are affine in $x$, and $\succeq_\mathcal{K}$ denotes the order induced by the convex cone $\mathcal K$. Moreover, the pairs $(x_\mathrm{lb}, g_\mathrm{lb})$ denote lower bounds, $(x_\mathrm{ub}, g_\mathrm{ub})$ denote upper bounds, and $(x_\mathrm{cb}, g_\mathrm{cb})$ denote conic bounds.

S = casos.sdpsol('S','solver',struct('x',x,'f',f,'g',g,'p',p),opts)

initializes the SDP solver named 'S' using the convex solver 'solver'. Options are provided as structure opts including optional fields opts.Kx and opts.Kc describing, respectively, the cones $\mathcal K_x$ and $\mathcal K_c$. See Convex cones for details.

sol = S('lbx',lbx,'ubx',ubx,'cbx',cbx,'lbg',lbg,'ubg',ubg,'cbg',cbg)

evaluates the SDP solver S providing (optional) arguments to describe $x_\mathrm{lb}$, $x_\mathrm{ub}$, $x_\mathrm{cb}$ and $g_\mathrm{lb}$, $g_\mathrm{ub}$, $g_\mathrm{cb}$.

Note

When using casos.sdpsol with Mosek, PSD decision variables must be constrained to the PSD cone on Kx, not Kc. This is because Mosek expects variables in triangulized (lower-triangular) form, while CaΣoS internally represents them in vectorized form. Symmetry is therefore enforced on Kx rather than on Kc. For further details, see Question #126.

Low-level interface

The low-level interface solves conic problems of the form

$$\begin{array}{l c r} \min & \frac{1}{2} x^\top H x + g^\top x, & x = (x_\mathrm{l}, x_\mathrm{c}) \\\ \text{s.t.} & a_\mathrm{lb} \leq A_\mathrm{l} \, x \leq a_\mathrm{ub}, & A_\mathrm{c} \, x \succeq_{\mathcal K_c} a_\mathrm{cb} \\\ \text{and} & x_\mathrm{lb} \leq x_\mathrm{l} \leq x_\mathrm{ub}, & x_\mathrm{c} \succeq_{\mathcal{K}_x} x_\mathrm{cb} \end{array}$$

where $\succeq_\mathcal{K}$ denotes the order induced by the convex cone $\mathcal K$. Moreover, the pairs $(x_\mathrm{lb}, a_\mathrm{lb})$ denote lower bounds, $(x_\mathrm{ub}, a_\mathrm{ub})$ denote upper bounds, and $(x_\mathrm{cb}, a_\mathrm{cb})$ denote conic bounds.

S = casos.conic('S','solver',struct('h',hs,'a',as),opts)

initializes the conic solver named 'S' using the convex solver 'solver', where hs and as are sparsity patterns for $H$ and $A = (A_\mathrm{l}, A_{\mathrm c})$. Options are provided as structure opts including optional fields opts.Kx and opts.Kc describing, respectively, the cones $\mathcal K_x$ and $\mathcal K_c$. See Convex cones for details.

sol = S('h',h,'g',g,'a',a,'lba',lba,'uba',uba,'cba',cba,'lbx',lbx,'ubx',ubx,'cbx',cbx)

evaluates the conic solver S providing (optional) arguments to describe $H$, $A$, and $g$ as well as $a_\mathrm{lb}$, $a_\mathrm{ub}$, $a_\mathrm{cb}$ and $x_\mathrm{lb}$, $x_\mathrm{ub}$, $x_\mathrm{cb}$.

Examples in CaΣoS

Semidefinite Program

Consider Example 6.7.1 from the Mosek documentation.

The problem involves semidefinite and conic constraints with a linear cost function, that is:

$$ \begin{aligned} &\text{minimize} \quad && \left\langle \begin{pmatrix} 2 & 1 & 0 \\ 1 & 2 & 1 \\ 0 & 1 & 2 \end{pmatrix}, \bar{X} \right\rangle + x_0 \\ &\text{subject to} \quad && \left\langle \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix}, \bar{X} \right\rangle + x_0 &= 1, \\ &&& \left\langle \begin{pmatrix} 1 & 1 & 1 \\ 1 & 1 & 1 \\ 1 & 1 & 1 \end{pmatrix}, \bar{X} \right\rangle + x_1 + x_2 &= \frac{1}{2}, \\ &&& x_0 \geq \sqrt{x_1^2 + x_2^2}, \\ &&& \bar{X} \succeq 0. \end{aligned} $$

In CaΣoS, the implementation of this problem is rather straightforward.

Decision Variables

We start by defining the symbolic decision variables as casadi.SX objects with

x_bar = casadi.SX.sym('x_bar', 3, 3);
x = casadi.SX.sym('x', 3, 1);

Note that $\bar{X}$ is a $3\times3$-matrix. Both are decision variables, meaning we define the sdp-struct with

sdp.x = [x; x_bar(:)];

Note

All problem decision variables need to be defined as one vector.

With the syntax

opts.Kx = struct('lin', [3], 'psd', [3]);

we define that there are linear bounds on $x_0, x_1, x_2$ and that $\bar{X} \in \mathbb{R}^{3\times3}$ is a positive semidefinite matrix. Since only $n \times n$ matrices can be positive semidefinite, we write ...'psd', 3... instead of ...'psd', 9....

Finally, we need to define the linear upper and lower bounds on $x_0,x_1,x_2$ with

lbx = -inf(3,1); % no bound -> set to -inf
ubx = inf(3,1); % no bound -> set to inf

Cost and constraints

We define the three matrices (cost function and first two constraints) with

A = [2 1 0; 1 2 1; 0 1 2];
B = eye(3);
C = ones(3, 3);

which now allows us to use standard Matlab syntax to define the cost function

sdp.f = trace(A*x_bar) + x(1);   % cost

and the three constraints with

sdp.g = [    % constraints
    trace(B*x_bar) + x(1);
    trace(C*x_bar) + x(2) + x(3);
    [x(1); x(2); x(3)];
]

where the actual dimension is size(sdp.g)= [5 1] since the Lorentz cone constraint considers three variables which we define with the K_-notation:

% cones
opts.Kc = struct('lin', [2], 'lor', [3]);

Since we could potentially have several Lorentz constraints, [3] can correctly be defined as a list (best practice). Finally, we add the upper and lower bounds for the first two linear constraints

lbg = [1; .5];
ubg = [1; .5];

to force the equality constraints given in the problem.

Now we can build the solver and obtain the solution with

% Use high-level interface
S = casos.sdpsol('S','mosek',sdp,opts);
sol = S('lbx',lbx,'ubx',ubx,'lbg',lbg,'ubg',ubg);

Complete code

The complete code is given by:

% https://docs.mosek.com/latest/toolbox/tutorial-sdo-shared.html#example-sdo1

% Decision variables
x_bar = casadi.SX.sym('x_bar', 3, 3); % matrix
x = casadi.SX.sym('x', 3, 1); % vector

sdp.x = [x; x_bar(:)]; % List of decision variables

lbx = -inf(3,1); % lower bound
ubx = inf(3,1); % upper bound


opts.Kx = struct('lin', [3], 'psd', [3]);

% Cost and constraints
A = [2 1 0; 1 2 1; 0 1 2];
B = eye(3);
C = ones(3, 3);

sdp.f = trace(A*x_bar) + x(1); % cost

% constraints
sdp.g = [
    trace(B*x_bar) + x(1);
    trace(C*x_bar) + x(2) + x(3);
    [x(1); x(2); x(3)];
]

% Conic constraints
opts.Kc = struct('lin', [2], 'lor', [3]);

% Linear upper and lower bound
lbg = [1; .5];
ubg = [1; .5];

% Use SDP interface
S = casos.sdpsol('S','mosek',sdp,opts);
sol = S('lbx',lbx,'ubx',ubx,'lbg',lbg,'ubg',ubg);