Robust optimization using fmincon in Matlab?

Question

I am trying to implement the following optimization (from this paper) in Matlab using fmincon:

$\min_\omega\omega'\Sigma\omega$ subject to $\min_Ur_p \geq r_0$ where $\Sigma$ is a positive definite matrix and $\omega$ is a vector of weights that sum to 1. Also: $r_p=\alpha'\omega$ and $U$ can be thought of as the sphere centred at $\alpha$ with radius equal to $|\chi|\alpha$ for $\chi$ between 0 and 1.

The authors show that:

$\min_Ur_p=|\alpha||\omega|[cos(\phi)-\chi]$ where $\phi$ is the angle betweent the two vectors $\alpha$ and $\omega$.

Any ideas how to implement this using fmincon?

For those interested, I am trying to implement the technique in this paper-specifically see pages 6-7.

Thanks in advance!

Answer 1

Based on the discussion in comments on your original question, it appears that your problem can be formulated in the form required by fmincon as

$ \min \omega^{T} \Sigma \omega $

subject to

$ -\alpha^{T}\omega + \chi \| \alpha \| \| \omega \| \leq -r_{0} $

$ \sum_{i=1}^{n} \omega_{i} = 1$

where $\alpha$, $\Sigma$, $\chi$, and $r_{0}$ are all known constants. You haven't said it, but I wouldn't be surprised if you also wanted to add the constraint $\omega \geq 0$. Note that I had to multiply the greater than or equal to constraint by negative one to turn it into a less than or equal to constraint.

To solve this problem using fmincon, you need to define a MATLAB function for objective and another MATLAB function for the nonlinear constraint. One slightly tricky issue here is how to get the parameters $\chi$, $\alpha$, and $\Sigma$ into these functions. This can be done using global variables (the old fashioned way) or using anonymous functions (the way that the Mathworks now recommends.) There are instructions on how to do this with an example in "help fmincon"

Answer 2

Just adding to the answers above, untangling the notation bit to put it in standard robust optimization format.

You have $r_p = \alpha + rU $ where the uncertainty $U$ is known to satisfy $U^TU \leq 1$ and your radius happen to be parametrized as $\kappa \left|\alpha\right|$.

After performing the minimization over the uncertainty, i.e, eliminating the for-all operator in the uncertain constraint, it boils down to a standard quadratic problem with a convex norm-constraint. Although fmincon most likely will work, it is not really the best tool for the task, since this is a very particular problem class for which there are dedicated extremly efficient solvers.

It could also be interesting to know that the problem easily generalizes to other norm-balls on the uncertainty.

For fun, here is some experimental code in the MATLAB Toolbox YALMIP

To begin with, solve the problem as manually derived above. If you have an efficient second-order cone programming solver installed, it will use that. Otherwise, it will probably use fmincon

n = 10;
alpha = randn(10,1);
S = randn(10);S = S'*S;
kappa = 0.01;
r0 = .01;
w = sdpvar(n,1);
Objective = w'*S*w;
Budget = [w >= 0, sum(w)==1];
Robust = [w'*alpha - kappa*norm(alpha)*norm(w) >= r0];
solvesdp([Budget,Robust], Objective)
double(w)

Actually, YALMIP has a built-in framework for deriving these kind of robust models automagically. Hence, the following problem derives the result in the paper and solves the problem

w = sdpvar(n,1);
U = sdpvar(n,1);
rp = alpha + kappa*norm(alpha)*U;
Objective = w'*S*w;
Budget = [w >= 0, sum(w)==1];
Uncertainty = [rp'*w >= r0, uncertain(U), U'*U <= 1];
solvesdp([Budget,Uncertainty], Objective)

This allows us to try alternative models, for instance a box-constrained uncertainty (this leads to a simple quadratic program, so fmincon is absolutely not the way to go)

w = sdpvar(n,1);
U = sdpvar(n,1);
rp = alpha + kappa*norm(alpha)*U;
Objective = w'*S*w;
Budget = [w >= 0, sum(w)==1];
Uncertainty = [rp'*w >= r0, uncertain(U), -1 <= U <= 1];
solvesdp([Budget,Uncertainty], Objective)

If you have an efficient mixed-integer SOCP solver installed (such as gurobi, cplex or mosek), you can for instance add combinatorial structures, such as a cardinality constraint on your positions

solvesdp([Budget,Uncertainty,nnz(w)<=3], Objective)

So, by looking at the problem a bit broader, I hope you realize there is a lot more fun that can be done, than simply throwing a general purpose nonlinear solver at a particular problem definition.

Answer 3

I'd recommend using different MATLAB routines instead, if you'd like to take the extra step of making your implementation faster, and are willing to install additional software.

Your problem looks like convex program that could be posed as a disciplined convex program, in which case using a solver like CVX has the potential to automatically transform the problem into a form more amenable to efficient solution (a semidefinite program, or a second-order cone program, perhaps).

That said, fmincon with the formulation suggested by Brian Borchers is probably the simplest and easiest implementation to get started. Your problem looks convex, and so it should not be too sensitive to choices of initial guesses.