LP file

I have a .lp file which takes very long time (more than 40 hours
without finish) when it is run on CPLEX's Interactive Optimizer
version 10.2.

As I only have CPLEX v.10.2, I was wondering if someone could run
the .lp file for me on the latest version of CPLEX (or Gurobi or other
solvers), and then tell me the results.

After you finish running the case, please tell me:
    (1) The software you use and the software version.
    (2) How long does it take to finish?
    (3) The spec. of the computer that are used to run the .lp file,
and whether the multi-thread function is turned on.
    (4) The final values of all the variables.
    (5) Why the .lp file takes so long?

This is a MIP model. An LP file (or MPS file) is a bad way to look at a MIP problem and give advice about better formulations. Difficult MIP models should always be implemented in a modeling language so you can look at the model, actually understand it and think about it.

Another good example how not to disseminate models is http://www.gamsworld.org/performance/index.htm. The models are totally obscured so they loose any meaning, and don’t allow any reasoning why a certain model gives performance problems for certain solvers.

When we want to learn something about the performance of a model I would argue that it is not helpful to use a format (LP/MPS) that hides the structure and that prohibits reformulations. This is especially the case for discrete models, but also for non-linear models.

Max Likelihood problems

1. I have this maximum likelihood problem. When I maximize the
likelihood function, I get some results; when I take the log of the
likelihood function and maximize it, I get different results. I
suppose this is due to the other problem I seem to have, which is too
many local maxima. Is this the typical case? Or is there something
fundamentally wrong in the way I am posing the problem?

2. Econometricians and statisticians like to use log-likelihood
functions instead of likelihood functions. The official argument is
that "it is easier to manipulate/optimize". Yes, it is easier to
derive, for humans. But it transforms a function that is bounded into
something unbounded. If the optimization is going to be made by a
computer does that really help or does that make it more difficult for
the computer/algorithm? I want to do whatever is the more reliable
technique. I hope this is not one of those "case-by-case" things.

Some points:

1. Taking the log can make the likelihood function better behaved. Some of the very large values that can be formed in sub-expressions can often be prevented this way.
2. Many good NLP solvers allow for bounds on the variables that allow you to steer the NLP solver away from the regions where the function can not be evaluated.
3. Use an alternative method to find a good starting point, e.g. a method of moments estimate.
4. Some of these problems can be difficult. Mixture models are a good example.

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)

GAMS: GDXXRW new option CheckDate

GDXXRW is a tool to exchange (read/write) data between GDX files and Excel spreadsheet files. The option CheckDate is helpful to speed up thing when importing data from Excel. If the GDX file is newer than the XLS file it will skip the actual data transfer. As this is COM based this step is somewhat expensive even fro small amounts of data (there is some overhead involved in starting Excel).

GAMS: log file problems

When running a large NLP in the GAMS/IDE with the new 23.1 version, I noticed that the log from KNITRO is buffered. This is not ideal: a log line should appear right away. Also the log from XA is completely disappeared since a few versions. In both cases the result is staring at a blank window when the job is running. The log is really the user-interface for a solver, so it would be good if a little bit more care is devoted to this.

Formulation: sum{i in I} binvar[i]*intvar[i]<=C

> How to linearize sum{i in I} binvar[i]*intvar[i] ≤ C where binvar is a binary variable and intvar is an integer variable.

In general we can linearize y[i] = binvar[i]*intvar[i] as:

0 ≤ y[i] ≤ intvar.up[i]
y[i] ≤ binvar[i]*intvar.up[i]
y[i] ≤ intvar[i]
y[i] ≥ intvar[i] - intvar.up[i]*(1-binvar[i])

In this case we can drop the ≤ inequalities, so we have:

0 ≤ y[i] ≤ intvar.up[i]
y[i] ≥ intvar[i] - intvar.up[i]*(1-binvar[i])
sum{i in I} y[i] ≤ C

Formulation c[i]=min(a[i],b[i])

> How to formulate c[i]=min(a[i],b[i]) for use with Cplex. A and b are integer.

There is not much we can exploit knowing that a, b are integer. In this formulation I use an extra binary variable bv[i] that implements the OR-functionality. The AMPL/Cplex formulation can look like:

var bv{I} binary;
e1{i in I}: c[i] <= a[i];
e2{i in I}: c[i] <= b[i];
e3{i in I}: bv[i]=0 ==> c[i]=a[i] else c[i]=b[i];

e3 implements c[i] = a[i] or c[i] = b[i] and e1,e2 make sure the minimum of the two is chosen.

GAMS: x**2 not allowed in QCP

> I have the the following model:

positive variable x;
free variable z;
equation obj;

obj.. z =e= x**2;
model m/all/;

solve m minimizing z using qcp;

> it is refused by GAMS.

I would consider that a bug. You can use any of the following alternative formulations:

z =e= x*x;
z =e= sqr(x);
z =e= power(x,2);

Formulation: c = Min(a,b) for binary a,b

> How to formulate c=min(a,b) where a,b are binary in a MIP.

This is like c = a×b. See: http://yetanothermathprogrammingconsultant.blogspot.com/2008/05/multiplication-of-binary-variables.html.

I.e.

c ≤ a
c ≤ b
c ≥ a + b − 1
c in [0,1] (continuous)
a,b in {0,1} (binary)

Similarly c = max(a,b) can be modeled via the constraints:

c ≥ a
c ≥ b
c ≤ a + b
c in [0,1] (continuous)
a,b in {0,1} (binary)

Formulation: Multiplication of integer variable with binary variable

> How to linearize z = x × b with x integer and b binary?

Use the formulation from http://yetanothermathprogrammingconsultant.blogspot.com/2008/05/multiplication-of-continuous-and-binary.html.

GAMS: Mystery Memory

Q: When I use profiling I sometimes see memory usage I cannot explain.

I suspect this is related to the following. Consider the model:

set i/i1*i500/;
alias (i,j,k);
parameter a(i,j,k);
a(i,j,k)$(uniform(0,1)<0.1) = uniform(0,1);

parameter b1(i),b2(j),b3(k);
b1(i) = sum((j,k),a(i,j,k));
b2(j) = sum((i,k),a(i,j,k));
b3(k) = sum((i,j),a(i,j,k));

When we run this with profiling on we see:

----      1 ExecInit                 0.000     0.000 SECS      3 Mb
----      4 Assignment a            20.779    20.779 SECS    410 Mb  12501427
----      7 Assignment b1            0.999    21.778 SECS    410 Mb    500
----      8 Assignment b2           11.544    33.322 SECS  1,216 Mb    500
----      9 Assignment b3           12.948    46.270 SECS  1,216 Mb    500

The assignment to a is somewhat expensive: we can do this faster (see http://yetanothermathprogrammingconsultant.blogspot.com/2009/07/fill-parameters-with-random-sparse-data.html). The memory use seems correct: 407 MB for 12501427 elements, which is about 34 bytes per nonzero element. We use the 64-bit version of GAMS. With the 32-bit we will see smaller numbers: 304 MB or about 25 bytes per element. (That explains why 64-bit versions are not faster as 32-bit version: they have to drag more Mb around as pointers use up more memory: 8 bytes vs. 4 bytes). This is all ok.

The assignment of b1 is also ok. It is fast and only adds 500 elements, so it does not even register in the total memory usage.

The assignment of b2 is strange. It is understandable it is slower: the summation is using indices in the “wrong” order. The sudden spike in memory usage is less obvious. The additional hefty 806 MB (approximately twice what is needed to store parameter a itself) is used internally to create faster access to parameter a. This automatic reordering option can be turned off with OPTION SYS11=1;. If we use this option the profile report looks like:

----      1 ExecInit                 0.000     0.000 SECS      3 Mb
----      1 Other                    0.000     0.000 SECS      4 Mb
----      5 Assignment a            20.467    20.467 SECS    410 Mb  12501427
----      8 Assignment b1            0.998    21.465 SECS    410 Mb    500
----      9 Assignment b2            7.925    29.390 SECS    410 Mb    500
----     10 Assignment b3           55.911    85.301 SECS    410 Mb    500

This report has the memory profile we expected. Note that the assignment b2 is actually faster without reordering. The assignment b3 however shows why reordering though expensive (in memory and CPU time) can really pay off. Note: the “Other” in the table is related to the handling of the OPTION SYS11.

So if you see unexpected large spikes in memory usage, this may well be related to the automatic reordering of symbols that can take place when assignments are not in the “correct” index order.

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.

Dense Input Format: does not work for large problems

In http://www.iging.com/optimization/ yet another input format for LP’s is proposed. I don’t think this particular format will work in practice: it stores all elements, even the zero ones. Large LP’s have typically a large number of nonzero elements. Lets take the medium sized LP indus89 (m=2,726, n= 6,570) from the GAMS model library:

Representation	Size (bytes)
GAMS: indus98.gms	262,612
MPS: fixed MPS generated by convert	1,634,841
proposed format (estim.)	143,359,504

For the large test model described in http://yetanothermathprogrammingconsultant.blogspot.com/2009/06/large-easy-lps-are-modeling-systems.html (m=1e6, n=1.2e6) we have:

Representation	Size (MB)
MPS: fixed MPS generated by convert	180
proposed format (estim.)	9,155,308

Anything dense will never work for large models (“dense does not make sense”). I doubt this will become a standard format. Besides the inability to store large problems, it seems that ASCII formats are more popular, not in the least because they are platform independent and because you can actually look at the content.

GAMS: Fill parameters with random sparse data

For testing I need sometimes very large, but very sparse data to play with. To generate random data, the easiest is something like:

a(i,j,k)$(uniform(0,1)<fracsparse) = uniform(-100,100);

However this will loop over every (i,j,k). If very sparse this can be become relatively expensive. It is faster in that case to loop over the expected number of nonzero elements, and place them randomly in the parameter. For this we use the lag operator for indices.

The two methods are not exactly the same: the second method can never generate more than nn nonzeroes, while the first method can. However for testing algorithms in GAMS the generated data is often “random enough” to give an indication about performance etc.

$ontext

   fill a(i,j,k) and b(i,j,k) with random sparse data.
   Both sparsity pattern as the data should be random.

   We demonstrate two methods. The second method is efficient if
   the parameters are *very* sparse.

$offtext

$set n 500
set i /i1*i%n%/;
scalar n /%n%/;
alias (i,j,k);


scalar fracsparse 'fraction of nonzero elements' / 0.0001 /;


*
* method 1
* this is slow because we loop over (i,j,k)
*

parameter a(i,j,k);
a(i,j,k)$(uniform(0,1)<fracsparse) = uniform(-100,100);


*
* method 2
* loop over number of nonzero elements only
*

scalar nn; nn = fracsparse*n*n*n;
display nn;

parameter b(i,j,k);
scalar kk;
loop(i$sameas(i,'i1'),
  for(kk=1 to nn,
     b(i+uniformint(0,n-1),i+uniformint(0,n-1),i+uniformint(0,n-1)) = uniform(-100,100);
  );
);


parameter statistics(*,*);
statistics('a','sum')   = sum((i,j,k),a(i,j,k));
statistics('a','count') = sum((i,j,k)$a(i,j,k),1);
statistics('b','sum')   = sum((i,j,k),b(i,j,k));
statistics('b','count') = sum((i,j,k)$b(i,j,k),1);
display statistics;

The timings are as follows:

	seconds
assignment a	18.237
assignment b (including loops)	0.031

Note1: the outer loop for b is not doing anything. It is looping over one element. The reason is we cannot use a lead or lag with a single element b(‘i1’+…).

Note2: This is not completely original: I think I have seen similar code in some model from Tom Rutherford. Sorry: I don’t have a reference for that.