MSF Solver Foundation API Question

Hi,
I am trying to port an Excel sheet that uses Solver the add-in from Frontline Systems to a C# program that uses the Solver Fondation API.
The sheet uses solver to find the combination of stocks in a portfolio for a given variance.
I have followed the Markowitz sample that comes with the express installation but one of the constraints that have to be added is that the portfolio's variance should be a certain value. The variance of the portfolio is calculated the following way:
SQRT(MMULT(MMULT(D38:L38;$B$10:$J$18);TRANSPOSE(D38:L38)))*SQRT(252)
where D38:L38 is the allocation vector (the result) and $B$10:$J$18 is the covariance matrix of the available stocks.
Could you please help me with how you would express this constraint in code?
thank,

I believe this ends up as a nonlinear constraint. In OML-like notation (I invent the Sqrt function here):

  Sqrt[Sum[{i,I},{j,I},x[i]*covar[i,j]*x[j]]] * Sqrt[252] == SomeValue

This shows an obvious advantage of using MSF/OML compared to Excel formulas: it is much easier to write readable constraints. Excel's Solver and its bigger brother from Frontline Systems can handle such nonlinear constraints through its GRG2 based NLP solver, but MSF does not support this. MSF supports convex QP's, but I don't think this can be reformulated into a linearly constrained QP.


Notes:

1. the lack of Sqrt is not a show-stopping problem, as we could write:
   
     Sum[{i,I},{j,I},x[i]*covar[i,j]*x[j]] == SomeOtherValue



2. For large problems this construct can be improved somewhat, first by using symmetry of Covar, and depending on the modeling system/solver by reformulating the quadratic form. As in your case the number of instruments is small, there is no need for that.



3. For very large problems sometimes a linear approximation is used (easy if minimizing variance, actually not so easy in this case).

MSF CSP: Tiling Squares (2)

In this post I displayed a CSP model for Microsoft Solver Foundation that solves the Tiling Squares problem. It was observed that a symmetry breaking constraint was really helpful. Basically it forced similar squares to be number in the x-axis direction. Of course a small extension would be to say: if same squares have the same x-coordinate, number them in the y direction. The symmetry breaking constraints can look like:

// this is to reduce some symmetry
Foreach[{i,T},{j,i+1,T},Implies[b[i]&b[j]&(Side[i]==Side[j]),x[j]>=x[i]]],
Foreach[{i,T},{j,i+1,T},Implies[
            b[i]&b[j]&(Side[i]==Side[j])&(x[j]==x[i]),
            y[j]>=y[i]]]

This reduces the calculation time further:

A bigger improvement can be achieved by observing that we can force lower numbered squares to be added first.

E.g. squares 9,10 and 11 have size 4. If we turn on square 10 we want to require that square 9 is also used. This can be implemented with the constraint:

FilteredForeach[{i,1,T},Side[i]==Side[i-1],Implies[b[i],b[i-1]]] 

This has much effect on the computation time:

 

If we add numbers to the squares we see how the numbering follows these rules:

The reference sequence A036444 shows a(19)=2. I.e. we can use just a set with a maximum of two identical tiles to cover an area of (19 × 19). Indeed we can verify this with this model. The data and results are:

MSF CSP: Tiling Squares

The tiling problem is not very amenable for solving by MIP solvers (see tiling.pdf). In this problem we try to cover an (n × n) square with smaller squares. In a MIP we need big-M formulations to express the non-overlap condition. With CSP we can write this down more directly and conveniently. Instead of using something like:

we can write using Microsoft Solver Foundation:

Foreach[{i,T},{j,i+1,T},
    x[i]>=x[j]+Side[j] |
    x[i]+Side[i]<=x[j] |
    y[i]>=y[j]+Side[j] |
    y[i]+Side[i]<=y[j]
]

Here the vertical bar means OR. A small (7 × 7) problem can be implemented quickly as follows:

Here we already used knowledge about which tiles to use. We only needed to determine where to place them. The real problem is slightly more complicated: we don’t know which tiles to use in advance. Below, we are looking for a configuration for the (9 x 9) problem where no tile can be used more than three times. This is number h(9) in Integer Square Tilings and a(9) in sequence A036444. To formulate this problem we add boolean variables b(i) indicating if tile i is being used.

Also we added an extra constraint to reduce symmetry (this caused significant speed-up). The solution times with and without the symmetry breaking constraint are:

Without symmetry breaking constraint	With symmetry breaking constraint

I did not make an exhaustive examination of all solver options, so there may be better performance possible. The complete model looks like:

Model[
    Parameters[Integers,N=9],   // size of square to fill
    Parameters[Integers,T=24],   // number of tiles
  
    Parameters[Sets,Tiles],   
    Parameters[Integers,Side[Tiles]],
       
    Decisions[Booleans,b[Tiles]],
    Decisions[Integers[0,8],x[Tiles],y[Tiles]],
   
    Constraints[
  
         Foreach[{i,Tiles},Implies[b[i],x[i]<=N-Side[i] & y[i]<=N-Side[i]]],
  
         // Note: Use T instead of Tiles
         Foreach[{i,T},{j,i+1,T},Implies[b[i]&b[j],
             x[i]>=x[j]+Side[j] |
             x[i]+Side[i]<=x[j] |
             y[i]>=y[j]+Side[j] |
             y[i]+Side[i]<=y[j]
         ]],
  
        Sum[{i,Tiles},Side[i]^2 * AsInt[b[i]]] == N^2,
          
        // this is to reduce some symmetry
        Foreach[{i,T},{j,i+1,T},Implies[b[i]&b[j]&(Side[i]==Side[j]),x[j]>=x[i]]]
    ]
  
]

We use the Implies[] construct here to implement if-then. Note that we used T instead of Tiles when we needed to operate on the upper-triangle of the (i,j) pairs. In this case we could not use Tiles as OML then assumes the indices are non-numeric (well, they are actually integers).

Update: more symmetry breaking constraints and their effect on solution times are discussed in this follow-up.

Fixed LP after MIP

Hi everybody,

After I solve a problem as an IP, I fix some of the variables
(complicating ones) to integer values of optimal solution to the IP
and solve the problem as a LP. The objective function only depends on
the fixed variables in both problems and there are constraints which
will be binding also those that will not. The problem is highly
degenerated.

The problem now is that I can not get any non-zero values for dual
variables using getDual or getDuals.

Anyone likes to put a comment on?

Cheers,

GAMS does this with all its MIP solvers to produce marginals for the rows and columns.

At the end of your fixed LP you should have a complete basis with m basics and n−m non-basics. The (non-binding) rows are likely to be basic, hence your zero's. The fixed integer variables are likely to be non-basic at bound, so you should see some nonzero marginals there.

I am sure your solver can produce a solution listing or a basis file so you can have a look where the non-basics end up and their marginal value (the dual value or reduced cost).

AMPL/Cplex: No solver needed

 

Hello,

Anyone can help?
I have run the following model using AMPL code with CPLEX 11.2 solver
in window xp. The program run well and take a couple seconds for the
small size of problem. But once the problem size is more than 200 (of
j and k), the solve time increases dramatically. When the problem size
of set j and set k is 900 ,  it took so long ( longer than 10 hrs) to
find the solution. I wonder, since the model looks simple with only 2
variables per constraint. I'm not sure did I write the AMPL code wrong
or there is something wrong about this. (However, I checked the
solution from CPLEX after it done. The solution is correct). Thank you
very much.

Best regards,
Alex

The Mathematical Model is:

Objective: Min Delta
S.t. (1) Xj *(Bj + Sum(i) M(ICij + TC1ij + TC2ijk)) <= Delta,     For all j, k

(2) Sum(i) Xj = 1

(3) Xj = {0, 1}

The AMPL code is:

set stag ;
set sup;
set w;

param tc1{i in sup,j in stag}>=0;
param b{j in stag}>=0;
param ic {i in sup,j in stag}>=0;
param M integer;
param tc2{i in sup,j in stag,k in w}>=0;

var x{j in stag} binary;
var Delta;

minimize cost: Delta;
subject to delt{k in w, j in stag}:  x[j]*(b[j]+ sum{i in sup} (M*(ic[i,j]+ tc1[i,j]+ tc2[i,j,k]) )) <= Delta;
subject to c: sum{j in stag} x[j] = 1;

Ten hours is a long time for 900 binary variables in this rather simple model. The problem has many rows as a result of the way the max is formulated. This construct where we minimize the max is sometimes difficult for MIP solvers. Note that some MIP models are just very difficult to solve (example). However in this case, I believe this problem can be solved extremely fast without even using Cplex. First we calculate a[j,k] such that we can simplify the problem to:

minimize cost: Delta;
subject to delt{k in w, j in stag}:  x[j]*a[j,k] <= Delta;
subject to c: sum{j in stag} x[j] = 1;

The first constraint can be interpreted as: if (x[j]=1) then Delta = max_k a[j,k]. The complete model can now be restated as:

param a{j in stag, k in w} := b[j] + sum{i in sup} (M*(ic[i,j]+ tc1[i,j]+ tc2[i,j,k]));
param d{j in stag} := max{k in w} a[j,k];
param dmin := min{j in stag} d[j];
var x{j in stag} binary := if (d[j]=dmin) then 1 else 0;
var Delta := dmin;
display x,Delta;

(I left as an exercise to make sure there is only one x[j]=1 in case multiple d[j]’s assume the value dmin). If you really want to call Cplex (especially after shelling out big bucks for it), the following should be very fast:

var x{j in stag} binary;
var Delta;

minimize cost: Delta;
subject to c: sum{j in stag} x[j] = 1;
subject to c2: sum{j in stag} x[j]*d[j] = dmin;
subject to c3: Delta = dmin;

option solver cplex;
solve;

This formulation automatically deals with the problem of multiple d[j]’s being equal to the minimum.

MSF QP problem

There is some issue with MS Solver Foundation’s QP. First I thought it was related to a reported problem about not reporting the constant term in objective values. But this must be something else:

gives:

The MSF people reported they know about the problem and it will be fixed for the next release. A workaround is to simplify the objective to obj->x*y-x-y+1. In that case the objective value is not reporting the constant term, but the levels of the variables x and y are correct. 

GAMS: how to get a list of parameters in a GDX file

> How can I get a list of the parameters in the farm.gdx? just the parameters, only the names

$call gdxdump farm.gdx symbols noheader | awk "{ if ($4==\"Par\") print $2 }" > farm_parameters.txt 

This is using a Unix tool awk, but this also works under windows as GAMS provides a copy of a number of such tools.

How to linearize this function

I'm trying to solve an optimization problem via linear or mixed integer programming model. But, the objective function includes the following expression f(x)=1/x + 1/x^2 Please let me know how to formulate this function as linear or mixed integer programming model

There are several possibilities:

1. Solve using a piecewise linear formulation (e.g. using SOS2 variables).
2. Solve using as SLP (Successive LP). In this case for simplicity you may want to try just linearizing (Taylor approximation) the objective each cycle.
3. Use an NLP solver (e.g. MINOS, SNOPT or CONOPT). You could solve first as LP (eg by dropping this term) and use that solution as starting point.

Sparse data in MSF OML

The modeling language of the Microsoft Solver Foundation does not support sparse data. I.e. if you have a parameter a[i,j] you need to have values for all a[i,j]’s and you cannot skip the zero’s. In some cases this is a little bit of problem. Adding all zero records may be a lot of work and may take a large amount of space in your spreadsheet. Sometimes it is possible to work around this problem by adding extra sets and columns. Probably the easiest way to explain this is with a small example. A network optimization problem is a good example where a sparse data structure can help. In this case we have a max flow model with a small but sparse network. In GAMS we can write:

$ontext
 
   max flow network example
 
   Data from example in
     Mitsuo Gen, Runwei Cheng, Lin Lin
     Network Models and Optimization: Multiobjective Genetic Algorithm Approach
     Springer, 2008
 
   Erwin Kalvelagen, Amsterdam Optimization, May 2008
 
$offtext
 
 
sets
   i 'nodes' /node1*node11/
   source(i)   /node1/
   sink(i)    /node11/
;
 
alias(i,j);
 
parameter capacity(i,j) /
   node1.node2   60
   node1.node3   60
   node1.node4   60
   node2.node3   30
   node2.node5   40
   node2.node6   30
   node3.node4   30
   node3.node6   50
   node3.node7   30
   node4.node7   40
   node5.node8   60
   node6.node5   20
   node6.node8   30
   node6.node9   40
   node6.node10  30
   node7.node6   20
   node7.node10  40
   node8.node9   30
   node8.node11  60
   node9.node10  30
   node9.node11  50
   node10.node11 50
/;
 
 
 
set arcs(i,j);
arcs(i,j)$capacity(i,j) = yes;
display arcs;
 
parameter rhs(i);
rhs(source) = -1;
rhs(sink) = 1;
 
variables
   x(i,j) 'flow along arcs'
   f      'total flow'
;
 
positive variables x;
x.up(i,j) = capacity(i,j);
 
equations
   flowbal(i)  'flow balance'
;
 
flowbal(i)..   sum(arcs(j,i), x(j,i)) - sum(arcs(i,j), x(i,j)) =e= f*rhs(i);
 
model m/flowbal/;
 
solve m maximizing f using lp;

The parameter capacity defines the network: where there is no capacity, we assume there no arc. We can calculate two-dimensional set arcs(i,j) to indicate which arcs exist. This can be used to simplify the flow-balance equations. Note that in GAMS only the variables that appear in the equations are actually part of the model. I.e. although we declare x(i,j), in the equations we only have x(arcs). This means we only have 22 variables x(i,j).

In OML we cannot use a sparse table for the capacities: we would need a complete table with 11x11=121 records. That would also mean that we have 121 variables. The trick is to have an explicit set arcs in the data table. As OML does not have 2-dimensional sets, I used parameters From[.], To[.] to describe the network. The complete data is below:

 
The OML model can now look like:


 
The node balance equation becomes a little bit unwieldy and less readable with this approach. This is a small, artificial example but note that quite a few models have some network component. It illustrates that leaving out sparse data handling and multidimensional sets to simplify a modeling language, although not a show stopper, comes with a cost in the form of additional complexity for the modeler. Note that in this model we still have a possibly large node table to maintain (the GAMS model handles this sparse).

AMPL vs. GAMS on radiotherapy treatment planning

 

hi all,
i would just like to solicit your advice regarding whether to get an
academic license for AMPL or GAMS. i'm currently doing a master's
thesis on cancer radiotherapy treatment planning (RTP), and i've yet
to decide between AMPL or GAMS for prototyping the models. some of my
decision criteria include language syntax, user interface, software
cost and available literature on RTP.
i think that the syntax for AMPL is better compared to GAMS because
the former more closely resembles the algebraic model formulations.
however, there seem to be more literature on RTP in GAMS and there's
also an integrated development environment for it. thanks for your
help.

 

Both systems are very good in what they are doing and widely used, so you cannot really go wrong with either choice. I would probably suggest to add extra points for the modeling system that is used by your colleagues and collaborators. That makes exchanging models and data easier and also is easier when discussing problems, tricks, issues etc.

Bob Fourer (AMPL) answered:

It's hard to find someone who can give equally expert advice on two competing
systems, as once you become familiar with one of them you don't usually have
much incentive to keep learning about the other.  But here are a few comments
from my hardly unbiased view.

AMPL was designed with the idea of being much closer to mathematical notation
and generally much more natural to use than GAMS, and it's superior on that
score.  A GAMS model typically relies on more special conventions and
reformulations than its AMPL counterpart; a case in point is the often extensive
use of the GAMS $ operator to impose various conditions.  Also, the IDE
notwithstanding, GAMS is fundamentally more of a batch system whereas AMPL
offers a more flexible option of interactively exploring models and results.
Finally while in certain areas GAMS is established through long use, still I see
modelers in these areas choosing AMPL, particularly when they are undertaking
new projects that do not depend on existing GAMS models.

In my opinion AMPL and GAMS are closer in practice than suggested here (e.g. where you use $ in GAMS, one would use : in AMPL). I actually slightly prefer the GAMS syntax when doing real work, as it is a little bit more compact and it is obvious where a summation ends (in AMPL this is based on operator priority, in GAMS a sum is visually bracketed by parentheses).