# Robust optimization using fmincon in Matlab?

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.

-
your statement of the constraint isn't clear at all. What do you mean by $\min_{U} r_{p}$ as a constraint? – Brian Borchers Feb 3 '13 at 1:32
It would really help if you defined all your terms: What are $\omega$, $\Sigma$, $U$, $r_p$, et cetera? Also, does it have to be fmincon? Based on the structure of the problem - which is not yet clear from the information you give - there could be better solvers available from Matlab. – Christian Clason Feb 3 '13 at 11:59
@Brian: I have corrected and made additions that should help to better explain the problem. – Geraldine Bailey Feb 3 '13 at 13:09
@Christian: I would like to use fmincon simply because that is the solver the authors of the paper use in their empirical analysis. However, I am open to other suggestions. – Geraldine Bailey Feb 3 '13 at 13:11
Is $\chi$ a fixed constant? – Brian Borchers Feb 3 '13 at 17:25

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"

-

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.

-
Hi Johan!!! I do appreciate this post as I implemented the optimization using fmincon and it takes quite some time to run! Close to 20 minutes. I will definitely try your suggestions...I have just downloaded the YALMIP toolbox. Thanks and kind regards. – Geraldine Bailey Feb 5 '13 at 21:53
If you have any further questions regarding YALMIP on this problem, post them at YALMIP groups – Johan Löfberg Feb 6 '13 at 10:15
Johan, we try to encourage people to stay on scicomp for a variety of reasons, including explicitly licensed CC BY content, customizable feeds, and a format designed to help people find answers to questions that have already been asked. I understand this approach doesn't work for everyone, but feel free to contact me if you have any questions about how we can better support YALMIP on here. – Aron Ahmadia Feb 6 '13 at 10:37
Aha, didn't know that. Thought it was the other way around, i.e., minimizing side discussion. – Johan Löfberg Feb 6 '13 at 11:17

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.

-