MSF Solver: DEA via Excel plug-in

he Excel plug-in allows only for one model. Data Envelopment Analysis models however consist of a number of small LP models: one for each DMU (Decision-Making Unit). However the individual LP models are very small. So it is reasonable to combine all small LP’s into one bigger LP. Basically we form:

max	+c1'x1	+c2'x2	...	+cK'xK
	A1x1=b1
		A2x2=b2
			...
				AKxK=bK

A basic CCR model can look like:

// DEA: combine all LPs into bigger LP
Model[
  Parameters[Sets,DMU,Inp,Outp],
  Parameters[Reals,x[DMU,Inp],y[DMU,Outp]],
  Decisions[Reals[0,Infinity],v[DMU,Inp],u[DMU,Outp]],
  Decisions[Reals,efficiency[DMU]],
  Constraints[
      // objectives:
      Foreach[{j0,DMU},efficiency[j0]==Sum[{o,Outp},u[j0,o]*y[j0,o]]],
      // normalize:
      Foreach[{j0,DMU},Sum[{in,Inp}, v[j0,in]*x[j0,in]] == 1],
      // limit:
      Foreach[{i,DMU},{j0,DMU},Sum[{o,Outp}, u[j0,o]*y[i,o]] <= Sum[{in,Inp}, v[j0,in]*x[i,in]]]
  ],
  Goals[Maximize[totaleff -> Sum[{j0,DMU},efficiency[j0]]]]
]

See also http://www.amsterdamoptimization.com/models/msf/oml.pdf. A similar trick was used in http://yetanothermathprogrammingconsultant.blogspot.com/2009/03/efficient-portfolio-frontier-in-msf.html.

Of course if you prefer to work at the API level there is no problem in solving a series of LP’s.

The Excel interface makes it easy to sort the results. In this case that is done with the Filter facility in Excel.

(click to enlarge)

MS Solver Foundation: Excel interface through COM/C# (2)

This is a follow-up on http://yetanothermathprogrammingconsultant.blogspot.com/2009/07/ms-solver-foundation-excel-interface.html. Here we explore an LP approximation to the QP models used in the original post. In the original model we used

risk = Σ w_t²/T

Here we will replace the quadratic expression by an absolute value:

risk ≈ Σ |w_t|/T

The big advantage of an LP model is that we can extend this to a MIP for more complicated modeling situations. Although some solvers support MIQP models, the MIQP algorithms are often not as fast as MIP algorithms. (Microsoft Solver Foundation does have a QP algorithm, but as it is based on an interior point algorithm it is not really suited to be used inside a branch-and-bound framework, so there is no standard MIQP capability).

The LP model is included in the DLL code as follows:

The absolute values are implemented using a standard variable splitting technique: w_t = y_t⁺ − y_t^- with y_t⁺,y_t^- ≥ 0. Then it follows that |w_t| = y_t⁺ + y_t^-.

Everything else can stay the same. The results are shown in the screen dump below. As you can see the LP offers a reasonable approximation: the portfolios for different trade-offs between risk and return are quite similar. (Click on the picture to enlarge).

This exercise shows that beyond using the MSF Excel plug-in, the .NET API’s in combination with the modeling language OML are a powerful tool to write applications and that these solutions can be easily integrated in an Office Environment (Access, Excel).

MS Solver Foundation: Excel interface through COM/C#

In preparation for a presentation on different modeling languages for a project at a financial firm, I implemented a demo QP/LP portfolio model with an Excel front-end and a .NET/C# DLL. 

MS Solver Foundation has a great tool to develop Excel based optimization models: the Excel plug-in.  One of the restrictions however is that only one single model can be handled this way. Sometimes we need to solve a series of models or different variants of a model. In this post I will show this can be achieved by creating a .NET DLL that uses the MS Solver Foundation API’s on one hand, and interfaces with Excel through a COM interface on the other hand.

Note that there are other possible approaches such as VSTO to handle this. VSTO allows you to create .Net applications inside an Office environment. The Microsoft Solver Foundation Excel plug-in is an example of a VSTO application. 

We used a portfolio problem as an example. Basically we have a multi-criteria problem:

minimize risk  x’Qx
maximize return  μ’x
subject to a budget constraint  Σx = 1

We attack this by solving a series of QP models

minimize λ x'Qx − μ'x
subject to a budget constraint  Σx = 1

for different values of λ.

The term for risk was rewritten as:

risk = Σ w_t²/T
w_t = Σ R_i,t x_i

where R are the mean adjusted returns. This separable formulation suggests an approximate LP formulation:

risk ≈ Σ |w_t|/T

The OML model for the QP formulation looks like:

Afterwards we can return x[i], Return and calculate the Risk term as (Obj+Return)*NumT/Lambda. We wanted to keep the quadratic term in the objective for obvious reasons.  Note that we need to substitute out {0} and {1} later so it holds values for the scalars Lambda and NumT. The substitution is done by String.Format, and we need to protect the other constructs like {t,T} by adding another pair of curly brackets: {{t,T}}.

The data for Mean and AdjRet are coming from an Access database. This is similar to http://yetanothermathprogrammingconsultant.blogspot.com/2009/05/ms-oml-model-using-c-and-data-from.html.

To generate a COM interface we use:

[ClassInterface(ClassInterfaceType.AutoDual)]

This will cause all public methods to become visible for VBA. For more information see: http://richnewman.wordpress.com/2007/04/15/a-beginner’s-guide-to-calling-a-net-library-from-excel/.

The QP model version is coded as:

This can be called from VBA as:

There is a little bit of extra code to retrieve the solution. At the end we collect solution data as:

In Excel it is easy to display this graphically:

The architecture of this demo app is as follows (click to enlarge):

The database not only stores the data, but also does some of the data manipulation: e.g. the mean adjusted returns are calculated in Access through queries.

CSP: indexing by a variable

A table lookup construct y=value[x] where x is an integer variable can be implemented in a MIP model as follows:

Parameters[Sets[Integers],I],
Parameters[Reals,value[I]],
Decisions[Integers[0,1],b[I]],
Constraints[
    Sum[{i,I},b[i]]==1,
    x == Sum[{i,I},i*b[i]],
    y == Sum[{i,I},value[i]*b[i]],
    ....

Here x and y are variables (they can be relaxed to be continuous as they are automatically integer – that is depending on value[] for y). Some experienced MIP users may recognize a SOS1 structure here. In a CSP model this can sometimes be expressed more succinctly as:

y == value[x]

This is also (and more coherently) discussed in the AMPL extensions document: http://www.ampl.com/NEW/FUTURE/logic.html#varsub. More info here: http://users.iems.northwestern.edu/~4er/WRITINGS/extmodcp.pdf. These references contain some motivating examples.

The Microsoft Solver Foundation group will add this to their OML language: http://code.msdn.microsoft.com/solverfoundation/Thread/View.aspx?ThreadId=1739. This will allow for very natural and intuitive formulations for some models.

MS OML model using C# and data from Access

Looking at http://code.msdn.microsoft.com/solverfoundation/Thread/List.aspx it is noticeable that many users have problems coding simple problems in C#. One possible reason is that they use C# to assemble the model. Although not extremely difficult, it is gives very unwieldy models: the signal-to-noise ratio in the code is small as you need lots of syntactic clutter just to specify all variables and constraints compared to a specialized modeling language. Large models can easily have dozens of blocks of variables and equations. In this post I want to emphasize an alternative that is somewhat underrated: it is possible to use OML directly in your C# application. That will immediately make the model much more compact and readable.

The second issue is that the data-binding is often not completely straightforward. Many questions are related to data binding. Below is the simplest solution I could come up with for the following architecture: the math programming model is a simple transportation model and all data is stored in an Access database. The goal is to provide a skeleton example that is both readable and simple while being useful as a starting point for larger, more complex applications. The advanced features of LINQ as used throughout the Solver Foundation documentation are largely geared towards SQL Server. Therefore I wanted to explore how a simpler database like Access could be handled. If Access is working, there is good reason to believe that any other major database will also work, as we are working with the lowest common denominator in some respects. Many databases are accessible through OleDb. In practice it may be a problem that all data has to come from the database: OML has no facilities for data manipulation. Large models often have large amount of data, which may need some form of processing (aggregation etc.). Even if your real database is say Oracle, it may be useful to use Access as front-end tool to perform these data manipulation steps.

The model is the trnsport.gms model from the GAMS model library. It is small and has a few small parameters. For more info see http://www.amsterdamoptimization.com/models/msf/oml.pdf. In OML the model looks like:

Model[
  Parameters[Sets,Plants,Markets],
  Parameters[Reals,Capacity[Plants],Demand[Markets],Cost[Plants,Markets]],

  Decisions[Reals[0,Infinity],x[Plants,Markets],TotalCost],

  Constraints[
     TotalCost == Sum[{i,Plants},{j,Markets},Cost[i,j]*x[i,j]],
     Foreach[{i,Plants}, Sum[{j,Markets},x[i,j]]<=Capacity[i]],
     Foreach[{j,Markets}, Sum[{i,Plants},x[i,j]]>=Demand[j]]
  ],

  Goals[Minimize[TotalCost]]
]

In C# we can do:

/// <summary>
/// Holds the OML model
/// </summary>
string strModel = @"Model[
      Parameters[Sets,Plants,Markets],
      Parameters[Reals,Capacity[Plants],Demand[Markets],Cost[Plants,Markets]],

      Decisions[Reals[0,Infinity],x[Plants,Markets],TotalCost],

      Constraints[
         TotalCost == Sum[{i,Plants},{j,Markets},Cost[i,j]*x[i,j]],
         Foreach[{i,Plants}, Sum[{j,Markets},x[i,j]]<=Capacity[i]],
         Foreach[{j,Markets}, Sum[{i,Plants},x[i,j]]>=Demand[j]]
      ],

      Goals[Minimize[TotalCost]]
   ]";

followed by:

SolverContext context;
context.LoadModel(FileFormat.OML, new StringReader(strModel));
Solution solution = context.Solve();
Console.Write("{0}", solution.GetReport());

This was easy and as short as can be. Now we need to get the data. The database is organized as:

We will use the tables Capacity and Demand and the Query Cost.  The data looks like:

To bind the data we use the following code:

/// <summary>
/// Solve the problem
/// </summary>
public void Solve()
{
    context.LoadModel(FileFormat.OML, new StringReader(strModel));

    foreach (Parameter p in context.CurrentModel.Parameters)
    {
        switch (p.Name)
        {
            case "Capacity":
                setBinding(p,"select plant,capacity from capacity",
                    "capacity", new string[]{"plant"});
                break;
            case "Demand":
                setBinding(p,"select market,demand from demand",
                    "demand", new string[]{"market"});
                break;
            case "Cost":
                setBinding(p,"select plant,market,cost from cost",
                    "cost", new string[]{"plant", "market"});
                break;
        }

    }

    Solution solution = context.Solve();
    Console.Write("{0}", solution.GetReport());

}

In each binding operation we specify:

1. The SFS parameter, which we retrieve from the CurrentModel
2. The query to be used against the database
3. The name of the data column
4. The names of the index columns (passed on as an array of strings)

The complete model looks like:

using System;
using System.Collections.Generic;
using System.Linq;
using System.Data;
using System.Data.OleDb;
using System.Data.Linq;
using System.Text;
using Microsoft.SolverFoundation.Services;
using System.IO;

namespace OML1
{
    class Trnsport
    {
        /// <summary>
        /// Called by the OS
        /// </summary>
        /// <param name="args"></param>
        static void Main(string[] args)
        {
            Trnsport t = new Trnsport();
            t.Solve();
        }

        /// <summary>
        /// Holds the OML model
        /// </summary>
        string strModel = @"Model[
              Parameters[Sets,Plants,Markets],
              Parameters[Reals,Capacity[Plants],Demand[Markets],Cost[Plants,Markets]],

              Decisions[Reals[0,Infinity],x[Plants,Markets],TotalCost],

              Constraints[
                 TotalCost == Sum[{i,Plants},{j,Markets},Cost[i,j]*x[i,j]],
                 Foreach[{i,Plants}, Sum[{j,Markets},x[i,j]]<=Capacity[i]],
                 Foreach[{j,Markets}, Sum[{i,Plants},x[i,j]]>=Demand[j]]
              ],

              Goals[Minimize[TotalCost]]
           ]";

        /// <summary>
        /// Connection string for MS Access 
        /// Use x86 architecture!
        /// </summary>
        string connection = @"Provider=Microsoft.ACE.OLEDB.12.0;Data Source=C:\projects\ms\OML1\OML1\trnsport.accdb;Persist Security Info=False;";

        /// <summary>
        ///  SFS
        /// </summary>
        SolverContext context;

        /// <summary>
        ///  Constructor
        /// </summary>
        public Trnsport()
        {
            context = SolverContext.GetContext();
        }

        /// <summary>
        /// get query result as DataSet
        /// </summary>
        /// <param name="connection">connection string</param>
        /// <param name="query">query as string</param>
        /// <returns></returns>
        private DataSet SelectOleDbSrvRows(string connection, string query)
        {
            DataSet ds = new DataSet();
            OleDbConnection conn = new OleDbConnection(connection);
            OleDbDataAdapter adapter = new OleDbDataAdapter();
            adapter.SelectCommand = new OleDbCommand(query, conn);
            adapter.Fill(ds);
            return ds;
        }

        /// <summary>
        /// Perform some magic to make sure the query output arrives in OML model.
        /// </summary>
        /// <param name="p">OML/SFS parameter</param>
        /// <param name="query">database query</param>
        /// <param name="valueColumn">column with values</param>
        /// <param name="IndexColumns">columns with indices</param>
        private void setBinding(Parameter p, string query, string valueColumn, string[] IndexColumns)
        {
             DataSet ds = SelectOleDbSrvRows(connection, query);
             DataTable dt = ds.Tables[0];
             p.SetBinding(dt.AsEnumerable(), valueColumn, IndexColumns);
        }

        /// <summary>
        /// Solve the problem
        /// </summary>
        public void Solve()
        {
            context.LoadModel(FileFormat.OML, new StringReader(strModel));

            foreach (Parameter p in context.CurrentModel.Parameters)
            {
                switch (p.Name)
                {
                    case "Capacity": 
                        setBinding(p,"select plant,capacity from capacity",
                            "capacity",new string[]{"plant"});
                        break;
                    case "Demand":
                        setBinding(p,"select market,demand from demand",
                            "demand", new string[]{"market"});
                        break;
                    case "Cost":
                        setBinding(p,"select plant,market,cost from cost",
                            "cost", new string[]{"plant", "market"});
                        break;
  
                }

            }

            Solution solution = context.Solve();
            Console.Write("{0}", solution.GetReport());

        }

    }
}

Some notes:

• This should work with almost any database. Just change the connection string accordingly.


• MS Access drivers are 32 bit so make sure you compile as 32 bit application. When targeting a 64 bit environment I got an exception about not being able to find an appropriate driver.


• The contents of the sets are derived from the parameter bindings: set elements are the union of the set elements used in the parameter binding.


• It may be useful to test the model beforehand using the Excel plug-in.


• Should use only one connection: make OleDbConnection conn a field of the object.


• How to write results back? This is often more complicated.

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.

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. 

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).

OML: MS vs Ketron

OML is the name of the modeling language (Optimization Modeling Language) of Microsoft Solver Foundation. I just remember seeing that TLA (Three Letter Acronym, itself a TLA) before. Ketron has OML defined as Optimization & Modeling Library. Ketron Management Science was famous for its MPSIII/Whizard mainframe based LP solver (competitor of IBM's MPSX). I believe Whizard used super-sparsity (this technology has largely disappeared together with the punch-card).

TSP in OML (MS Solver Foundation)

In theory it is easy to formulate a TSP when the all-different constraint is present. Let x[i] be the i^th city visited, an OML-like formulation could look like:

Model[

  Parameters[Integers,N=14],

  Parameters[Sets[Integers],city],

  Parameters[Reals,dist[city,city]],

  Decisions[Integers[0,N-1],x[city]],

  Constraints[

     Unequal[Foreach[{i,city},x[i]]]

  ],

  Goals[

    Minimize[Sum[{i,city},dist[x[i],x[i++1]]]]

  ]

]

This does not work completely. First, bounds can not contain a symbolic constant. Yuck: we need to write [0,13] as bounds:  Decisions[Integers[0,13],x[city]]. The objective has  a few more difficulties. There is no circular lead/lag operator in OML (like the -- or ++ operator in GAMS).  This can be solved by introducing a parameter next[city]. This can be calculated in Excel and imported in the OML model. Furthermore, we cannot use a variable as an index, so dist[x[i],x[next[i]]] is not a valid construct (some languages geared towards solving CSP problems allow this; it adds concise expressiveness to the language that may help the modeler). The error message indicates this is not related to OML per se but rather to solver support. Also it is not allowed to use something like FilteredSum[{j,city},j==x[i],...] (a condition in a FilteredSum cannot depend on a decision variable). The only thing I could think of was:

Minimize[Sum[{i,city},{j,city},{k,city},AsInt[j==x[i]]*AsInt[k==x[next[i]]]*dist[j,k]]]

Note that I used integers as set elements. I started with using data-binding through tables, using set elements {'city0','city1',...,'city13'}. This made the model somewhat more complicated, as x[i] is an integer. So the CSP formulation uses the simpler set {0,1,...,13} allowing us to operate on indices. (Some consider this bad modeling: in GAMS arithmetic on sets is actually discouraged, and needs a function ord()to convert the element — always a string in GAMS — to an integer). Unfortunately I was not able to solve the small 14 city example with the CSP solver using this formulation. (Even after adding City[0]==0 as constraint). To check the input, I also tried a MIP formulation with Gurobi. That solved this small instance just fine.

In the MIP formulation we needed to formulate:

FilteredSum[{i,city},i != 'city0', x['city0',i]]==1

From what I can see this is not possible: a set element not being an integer can not be used in OML. As a workaround I created sets by binding to some dummy parameters:

Sum[{i0,city0},{i,city2},x[i0,i]]==1

where city0 is a set containing only a single element 'city0', and city2 is a set {'city1',...,'city13'}.  The constraint

FilteredSum[{i,city},{j,city},i != j,x[i,j]] <= N

did not work as i and j are not numeric. A workaround could be to have a parameter num[city] which can be calculated in Excel and then imported. Then the condition can read: num[i]!=num[j]. I just used dist[i,j]>0 as condition, as that also can be used to exclude the diagonal (the Excel spreadsheet makes sure all diagonal distances dist[i,i] are zero). A more general approach would be to be able to introduce sparse 2-d sets (like one could do in AMPL and GAMS), but OML has only one dimensional sets as far as I know.

This exercise was really meant to explore the expressiveness of OML. The little model actually shows some interesting issues with OML and some possible workarounds. I do not want to suggest this is a good way to solve TSP's. Obviously these approaches are not suited for large problems. A cutting plane algorithm often is quite effective for slightly larger problems (see this 42 city problem). The excel plugin does not allow us to to implement such an algorithm: an OML model can only contain a single model and has no looping facilities. For the real large problems, you need to look elsewhere, such as Concorde.

MS Solver Foundation /Gurobi (2)

The MIP capabilities of the standard editions have been significantly reduced. In the previous version this was 2k variables and equations and now only 500 variables and equations. This is looks like a Gurobi limit:

Error:
The solver(s) threw an exception while solving the model - Gurobi throttle exceeded.
Visit www.gurobi.com/html/products.html to purchase an unlock license (HostID: 8AA3C6D2).


The obj seems to be counted as an equation:

Model[

  // generate simple MIP model with 500 variables and 500 constraints
  Parameters[Integers,N=500],
  Decisions[Integers[0,1],Foreach[{i,N},x[i]]],
    
 Constraints[ Foreach[{i,N},x[i] <= 1]  ],
 
  Goals[
    Maximize[Sum[{i,N},x[i]]]
    //Maximize[obj->Sum[{i,N},x[i]]] use to prevent funny name for goal
  ]

]

does trigger the above message. From the Excel plug-in I don't see an option to use the MS MIP solver. I hope I am wrong, but I suspect there is a class of MIP models that could be solved with the previous version of the Standard Edition that can no longer be solved with the current 1.1 release.

The workaround seems obvious: purchase the enterprise edition or purchase Gurobi.

Note that the limit is on any variable, i.e. 1 integer variable + 500 continuous variables makes the model too large.


Update: If your problem is in between these limits this can help. There is a way to point the excel plugin back to the MS MIP solver through a configuration file. For more info please contact the MSF people at http://code.msdn.microsoft.com/solverfoundation. They are very helpful and responsive.

MS Solver Foundation / Gurobi

MS Solver Foundation 1.1 is out. It includes Gurobi! (it has even become the default MIP solver). See http://code.msdn.microsoft.com/solverfoundation.

Indeed we see when running a MIP from the Excel plugin:

Social Golfer Problem in OML

I am recently involved in a practical application of a scheduling problem related to the Social Golfer Problem. The particular application is somewhat more complicated, but we could use the GAMS/Cplex model described here as a starting point.

For smaller instances of the pure Social Golfer Problem, it is convenient to use a CSP approach. Here is a formulation in OML (Microsoft Solver Foundation):

Model[
 
// N : number of golfers
// NG : number of groups
// T : number of rounds
// GS : group size (N/NG)
 
 Parameters[Integers,N=16,NG=4,T=5,GS=4],
 // Would prefer: GS=N/NG but OML does not allow (constant) expressions
 // in parameter statements
 
 Decisions[
  Integers[0,1],
   Foreach[{i,N},{g,NG},{t,T},x[i,g,t]]
 ],
 
  
 Constraints[
     // each golfer has to play each round
     Foreach[{i,N},{t,T},Sum[{g,NG},x[i,g,t]] == 1],
 
     // form groups
     Foreach[{g,NG},{t,T},Sum[{i,N},x[i,g,t]] == GS],
 
     // golfer i and j meet at most one time
     // Foreach[{i,N}, FilteredForeach[{j,N},i!=j, Sum[{g,NG},{t,T},x[i,g,t]*x[j,g,t]] <= 1]],
     // We exploit symmetry here:
     Foreach[{i,N}, Foreach[{j,i+1,N}, Sum[{g,NG},{t,T},x[i,g,t]*x[j,g,t]] <= 1]],
 
 
     // fix first round
     Foreach[{g,NG},{k,GS},x[GS*g+k,g,0]==1],
 
     // fix first golfer
     Foreach[{t,1,T},x[0,0,t]==1]
 ]
 
]

The condition on number of times players can meet can be implemented quite straightforwardly. In the MIP model we needed to linearize this constraint, but in the CSP framework we can use the nonlinear formulation directly. A big advantage of CSP over MIP is that many constructs can be modeling in a more natural, straightforward fashion. Note that we added some exploitation of symmetry: once we have checked how many times golfer i meets golfer j we no longer need to check golfer j meeting golfer i.

The last two constraints fix the first round and first player. This (without loss of generality) reduces the size of the problem.

This instance solved easily:

Solver Foundation Report       
Report Overview       
Date Time    =    1/18/09 10:48 AM
Model Name    =    golf.xlsx
Solver Selected    =    CSP
Total Time (ms)    =    677
Solver Completion Status    =    Feasible
       
Solver Execution Details       
Algorithm    =    Tree Search
Variable Selection    =    Dynamic Weighting
Value Selection    =    Forward Order

Note that we changed the default solver algorithm settings to speed up the solution times.

The most interesting instance with 32 golfers, foursomes and T=10 is too large for my version to solve:

Parameters[Integers,N=32,NG=8,T=10,GS=4],
 
 
Error:
Model size limit has been exceeded for this version of the product. Please contact Microsoft Corporation for licensing options.

Mona Lisa assignment in GAMS

This is the assignment problem described in Donald E. Knuth, "The Stanford GraphBase: A Platform for Combinatorial Computing", ACM Press, 1993. The gray scale data was generated by lisa(360,250,255,0,360,0,250,0,22950000) and is used to color (small) Excel cells.

The GAMS model is as follows:

$ontext
  
   Assignment problem for lisa(360,250,255,0,360,0,250,0,22950000)
  
   Erwin Kalvelagen, March 2002
  
   Reference: Donald E. Knuth, "The Stanford GraphBase:
        A Platform for Combinatorial Computing", ACM Press, 1993.
  
  
$offtext
  
  
set i 'rows' /0*359/;
set j 'columns' /0*249/;
  
$if exist input.gdx $call 'rm -f input.gdx'
$if exist output.xls $call 'rm -f output.xls'
 
 
$call '=gdxxrw.exe i=input.xls par=c rng=Data!A1 rdim=1 cdim=1'
 
parameter c(i,j);
$gdxin input.gdx
$load c
  
equations
    objective 'maximize highlighted parts'
    column(j) 'column sum'
    row(i)    'row sum'
;
 
variables
    z         'objective variable'
    x(i,j)    'pixels'
;
 
binary variable x;
 
objective..      z =e= sum((i,j),c(i,j)*x(i,j));
row(i)..         sum(j, x(i,j)) =l= 1;
column(j)..      sum(i, x(i,j)) =e= 1;
 
model m /all/;
solve m using rmip maximizing z;
 
set solution(i,j);
solution(i,j)$(x.l(i,j)>0.5) = yes;
execute_unload 'output.gdx',solution;
 
execute '=gdxxrw.exe i=output.gdx o=output.xls set=solution rng=Solution!A1 rdim=2 cdim=0'

This model is stored inside the Excel spreadsheet. The Excel front-end will export the model and the input data to a temp directory. Then GAMS is called, and finally the resulting solution data is read. To make it possible to use gdxxrw inside the GAMS model, we use the following algorithm:

1. Write GAMS model to %TEMP%\MonaLisa.gms
2. Write Data sheet to %TEMP%\input.xls
3. Call gams to execute %TEMP%\MonaLisa.gms
4. Copy %TEMP%\output.xls to Solution sheet
5. Draw points in Solution sheet with different color

Note that GDXXRW cannot read/write directly to the spreadsheet that has VBA code running, neither can we make it shareable. This algorithm bypasses this problem by copying input and output data to external workbooks.

The model is a large but easy RMIP: 600 equations and 90,000 variables. The free LP solver CBC can solve this quickly: 0.6 seconds. Note: the model is too large to solve with MS Solver Foundation Express and Standard Editions (see: note). Using the Check button in the MS Solver Foundation Excel plug-in seems to indicate that the Excel binding is somewhat slow. In addition the plug-in does not provide a progress window (see second picture below) which is very useful for long running MIP models. Although the GAMS/EXCEL/VBA route requires more programming, for some classes of models this approach seems more flexible than the MS Solver Foundation Excel plug-in solution. Of course for ultimate flexibility at the cost of even more programming effort you can use the MS Solver Foundation .NET API's.

MSF CSP: Domain Over Weighted Degree

The setting Domain Over Weighted Degree seems to be useful for a number of CSP models in the Microsoft Solver Foundation. 

Examples:

1. Magic squares example.
   
2. User model.

MSFCli.exe

Some issues I raised wrt MSFCLi.exe will be fixed in the next release of MS Solver Foundation.

Redeclarations in GAMS

> Erwin:
> To my surprise GAMS allows
> variable x,x;
> Any reason for this?

GAMS allows redeclarations of symbols unless there is a contradiction. I.e. one often sees:

variable x(i),y,z;

positive variable x;

This is perfectly legal GAMS. Not allowed is:

positive variable x;

negative variable x;

Your example is allowed as there is no contradiction. 

MS Solver Foundation: Excel/OML issue

I was sometimes loosing the model information when exiting Excel. Turns out that when you save the spreadsheet after closing it (Excel will prompt you if you want to save your changes) the OML model gets lost. (Somewhat a warped way of saving the workbook, but that is how my brain is wired.) If you manually save the model before closing Excel, the model will be saved without a problem. Microsoft could reproduce it and it is quickly fixed in the next version. (See issue report).

I had hoped to solve the assignment problem suggested by Knuth. It is supposed to find "interesting" spots on the Mona Lisa. Unfortunately, the standard edition of MS Solver Foundation is "throttled" too tightly to be able to solve this large but easy model. I am sure this is not a good business case to purchase the Enterprise version! (See limits).