I am a full-time consultant and provide services related to the design, implementation and deployment of mathematical programming, optimization and data-science applications. I also teach courses and workshops. Usually I cannot blog about projects I am doing, but there are many technical notes I'd like to share. Not in the least so I have an easy way to search and find them again myself. You can reach me at erwin@amsterdamoptimization.com.

## Tuesday, September 30, 2008

### Example MSFCli.exe

C:\projects\test>"c:\Program Files (x86)\Microsoft Solver Foundation\1.0.6.4960\Bin\MSFCli.exe" +type mps +verbose +log indus89.mps

===== Processing .\indus89.mps =====

Algorithm: Primal

CostingExact: SteepestEdge

CostingDouble: SteepestEdge

Variables: 6570 -> 6366 + 2363

Rows: 2727 -> 2564

NonZeros: 39490

Eliminated Slack Variables: 201

Phase 1 double Pivots: 4639

Phase 2 double Pivots: 1864

double Factorings: 138

Rational Factorings: 1

Degenerate pivots: 1367 (21.02 %)

{

x1 -> 114873.65553186,

x2 -> 1522.85610775999,

x3 -> 3460.00034188404,

x4 -> 479.152767198711,

x5 -> 22532.3511432455,

x6 -> 2490.55460455224,

....

x6530 -> 1,

x6550 -> 1,

x6570 -> 1

}

Zero variables: 3949

Positive variables: 2621

Negative variables: 0

Objective value: 'obj' Min Optimal -114873.65553186

===== Finished .\indus89.mps: 00:00:05.8776890 Optimal =====

The solver seems to feature some kind of rational arithmetic.

Not sure if this is related, but the Microsoft solver finds a correct solution to the Neumaier model http://www.gams.com/modlib/libhtml/badmip.htm.

### Microsoft Solver Foundation

- Simplex solver (LP,MIP)
- Interior Point solver (LP,QP)
- Constraint solver (CSP)
- Unconstrained non-linear solver (NLP)

Can be used as library via any .NET language, or input can be specified as MPS and a simple equation oriented format called OML. In addition, Excel bindings are available.

http://code.msdn.microsoft.com/solverfoundation

After downloading, documentation and samples are available in My Documents\Microsoft Solver Foundation.

## Monday, September 29, 2008

### GBD feasibility cut

*>Dear Mr. Kalvelagen,*

>

>I am trying to do a comparison between the performance

>of GBD and OA in the case of convex and nonconvex minlp

>problems.

>

>I am using the GBD-code from the script for MINLP methods

>that you have written and posted on your page, and comparing

>it with DICOPT-OA. My application is on transportation

>network design problem. The theory as highlighted in Floudas

>1995 about the difference between the two methods seems to

>match for the most part in my problem. However, I am facing

>infeasibilities in the sub-problem.

>

>Your GBD code on the portfolio optimization generates support

>functions only from the feasible primal. Can you give me an

>idea (if that is not much of a hassle for you) how a

>feasibility primal subproblem can be implemented in your

>GBD code to generate additional supports (from the infeasible

>subproblems), or if you have provided an example of such a

>case somewhere else I’d be grateful if you can direct me to that?

>

>I do find your webpage and the relevant files very useful when

>testing theoretical concepts through implementation in GAMS.

>

>Many thanks in advance for your time!

>

>I am trying to do a comparison between the performance

>of GBD and OA in the case of convex and nonconvex minlp

>problems.

>

>I am using the GBD-code from the script for MINLP methods

>that you have written and posted on your page, and comparing

>it with DICOPT-OA. My application is on transportation

>network design problem. The theory as highlighted in Floudas

>1995 about the difference between the two methods seems to

>match for the most part in my problem. However, I am facing

>infeasibilities in the sub-problem.

>

>Your GBD code on the portfolio optimization generates support

>functions only from the feasible primal. Can you give me an

>idea (if that is not much of a hassle for you) how a

>feasibility primal subproblem can be implemented in your

>GBD code to generate additional supports (from the infeasible

>subproblems), or if you have provided an example of such a

>case somewhere else I’d be grateful if you can direct me to that?

>

>I do find your webpage and the relevant files very useful when

>testing theoretical concepts through implementation in GAMS.

>

>Many thanks in advance for your time!

You can add feasibility cuts by forming a feasibility subproblem.

I do something like that for the dual subproblem in http://amsterdamoptimization.com/pdf/benders.pdf

A straightforward discussion can be found in Ullmann's Modeling and Simulation By Fritz Ullmann page 113, and 116.

Of course it may be possibly to bypass the problem by formulating an elastic version such that the subproblems are always feasible.

### Conditional generation

*> OK thanx for this information.*

> Did you have time to look at my "corrected" constraints? I'm really sorry for

> having written an error ...

> The correct conditions are now written:

>------------------------------------------

> if s=1 then t > b

> if s=-1 then t < -b

> Did you have time to look at my "corrected" constraints? I'm really sorry for

> having written an error ...

> The correct conditions are now written:

>------------------------------------------

> if s=1 then t > b

> if s=-1 then t < -b

*> if s=0 then -b <= t <= b*

*> where:*

> variables: t real in [-1,1] b real in [0,1]

> parameter: s integer in {-1,0,1}

> ------------------------------------------

>I just can't see how to write these conditions as linear constraints in AMPL ...

>Thanx again for your help or any advice which could guide me towards the solution ...

> variables: t real in [-1,1] b real in [0,1]

> parameter: s integer in {-1,0,1}

> ------------------------------------------

>I just can't see how to write these conditions as linear constraints in AMPL ...

>Thanx again for your help or any advice which could guide me towards the solution ...

param s = 1;

var t >= -1, <= 1;

var b >= 0, <= 1;

minimize obj: b; subject to

e1{1..1:s=1}: t >= b;

e2{1..1:s=-1}: t <= -b;

e3{1..1:s=0}: -b <= t;

e4{1..1:s=0}: t <= b;

Conditional generation of equations is an important feature in formulating real-world models. You may want to consider hiring a consultant with experience in AMPL modeling. That would speed up the implementation considerably.

## Wednesday, September 24, 2008

### GAMS Command Prompt

## Thursday, September 18, 2008

### $load x=x.l not allowed (GAMS bug)

execute_load 'merged.gdx',x=x.l;

$gdxin merged.gdx

$load x=x.l

GDXIN C:\projects\china\input\merged.gdx

5925 $load x=x.l

**** $490

**** LINE 19 INPUT C:\projects\china\input\report.gms

**** 490 Cannot load/unload this suffix

## Tuesday, September 16, 2008

### On solving a large system of nonlinear equation

_{i}

_{i}+w

_{i}

## Monday, September 15, 2008

### sub models

However this does not work well with the IDE. For instance the log of a submodel if called by

$call =gams submodel.gms

will be invisible and lost. There are few thing we can do to improve the situation. The settings I use are as follows:

`* flags to run submodels better under the IDE`

$setglobal ide "ide=%gams.ide% lo=%gams.lo% errorlog=%gams.errorlog% errmsg=1"

$call =gams submodel.gms %ide%

$if errorlevel 1 $abort submodel.gms has issues

This will solve most problems of calling gams via $call when using the IDE, and also works good when not running under the IDE.

The one remaining and somewhat annoying problem (bug) is that clicking on black lines of the submodel log, will cause the IDE to go to listing file of the main model. I.e. the wrong listing file is loaded. Clicking on the red lines work fine: the IDE will go the correct source file.

I think it would be an improvement if GAMS has better built-in support for executing submodels (compile time and execution time) where we don't have to worry about all these parameters and with better interaction with the IDE.

## Friday, September 12, 2008

### Restarts

*> I am dealing with a class of large scale LPs. The class has nice*

> structure though, let's say it has 1000 constraints. The objective

> function and the first 999 constraints are identical among all

> instances of the problems in this class. So everytime that user wants

> to solve a new problem it is only different from other problems in the

> 1000th constraint.

>

> I was wondering if this feature can be exploited to improve the speed

> of solving problems in this class. Is it possible that something be

> pre-computed and stored about the problem from the fixed part

> (objective and 999 constraints) and then everytime that user inserts

> the 1000'th constraint, that knowledge be used efficiently?

> structure though, let's say it has 1000 constraints. The objective

> function and the first 999 constraints are identical among all

> instances of the problems in this class. So everytime that user wants

> to solve a new problem it is only different from other problems in the

> 1000th constraint.

>

> I was wondering if this feature can be exploited to improve the speed

> of solving problems in this class. Is it possible that something be

> pre-computed and stored about the problem from the fixed part

> (objective and 999 constraints) and then everytime that user inserts

> the 1000'th constraint, that knowledge be used efficiently?

Absolutely. 1000 constraints is very small, but in general simplex based LP solvers have excellent restart facilities using the concept of a basis. Using the GAMS modeling language this even completely automated. I.e. the fragment:

`set scenario /s1*s5/;`

parameter growthq_s(scenario) /

s1 2

s2 2.2

s3 2.3

s4 2.5

s5 2.7

/;

parameter report(scenario,*);

loop(scenario,

growthq = growthq_s(scenario);

solve wsisn maximizing cps using lp;

report(scenario,'time') = wsisn.resusd;

report(scenario,'iterations') = wsisn.iterusd;

);

display report;

solves a series of related LP's: one coefficient in the model

changes here. Iteration count and cpu time are recorded and

they are displayed as:

` ---- 3642 PARAMETER report`

time iterations

s1 0.630 3704.000

s2 0.108 76.000

s3 0.070 76.000

s4 0.066 76.000

s5 0.069 76.000

so indeed the first solve is the most expensive one.

### Bug in gams

1 SET NUM /1*2/;

2 VARIABLE X(NUM);

3 X.L(NUM) = ORD(NUM);

4

5 X.FX('2') $(X.LO('2') EQ X.UP('2')) = 234;

6 DISPLAY X.LO,X.L,X.UP;

COMPILATION TIME = 0.000 SECONDS 3 Mb WEX228-228 Jul 26, 2008

---- 6 VARIABLE X.Lo

2 234.000

---- 6 VARIABLE X.L

1 1.000, 2 234.000

---- 6 VARIABLE X.Up

2 234.000

### GAMS/Xpress documentation bug

## Thursday, September 4, 2008

### Smooth approximations

- ifthen(x1>x2, y, z) = z + (y-z)*0.5*(1+tanh(k*(x1-x2)))
- ifthen(x1<x2, y, z) = y + (z-y)*0.5*(1+tanh(k*(x1-x2)))
- max(x,y) = ifthen(x>y,x,y) = y + (x-y)*0.5*(1+tanh(k*(x-y)))
- min(x,y) = ifthen(x<y,x,y) = x + (y-x)*0.5*(1+tanh(k*(x-y)))

## Wednesday, September 3, 2008

### Bug in ifthen (GAMS)

` 1 variable x;`

2 equation e1,e2;

3

4 scalar s /1/;

5

6 e1.. x$(s>0) =e= 0;

7 e2.. ifthen(s>0,x,0) =e= 0;

8

9 model m1 /e1/;

10 solve m1 using cns;

11

12 model m2 /e2/;

13 solve m2 using cns;

The models are the same. The ifthen is exogenous but GAMS does not recognize this. As a result the fist model solves fine and the second model leads to numeric problems:

** Error in Square System: Pivot too small.

Turns out the gradients are not correct for ifthen:

`---- e2 =E= `

e2.. (0)*x =E= 0 ; (LHS = 0)

Advice: don't use ifthen in equations as the derivatives of this function are wrong. If possible use $ conditions (exogenous if) or a smooth approximation (endogenous if). Smooth approximations for ifthen are discussed in this blog.

### Arithmetic on indices

*>Is it possible to write an objective function as*

>"obj.. PROFIT =e= sum((l,m,dig) $(ord(l) le p and ord(m) le p and ord(dig) le digmax ), x('(1+p*(ord(l)-1))','(1+p*(ord(m)-1))',dig)*ord(dig)); "

>if I only want to add particular x's to my objective function (indexed as 1+p*(ord(l)-1) and 1+p*(ord(m)-1) where p is a constant )?

>"obj.. PROFIT =e= sum((l,m,dig) $(ord(l) le p and ord(m) le p and ord(dig) le digmax ), x('(1+p*(ord(l)-1))','(1+p*(ord(m)-1))',dig)*ord(dig)); "

>if I only want to add particular x's to my objective function (indexed as 1+p*(ord(l)-1) and 1+p*(ord(m)-1) where p is a constant )?

*>In this case, gams gives me domain violation error for first and second index and i couldn't find a way to solve it*

*>(the domain is not violated mathematically for example for p=3 and digmax= 9).*

*>Is there a way to do arithmetic operations in indexes?*

Indices are always strings in GAMS. But you can do what you want as follows:

`alias (l,i), (m,j); `

obj.. PROFIT =e=

sum((l,m,dig,i,j) $(ord(l) le p and ord(m) le p and ord(dig) le digmax and ord(i)=(1+p*(ord(l)-1)) and ord(j)=(1+p*(ord(m)-1))),

x(i,j,dig)*ord(dig));

This is somewhat messy, but we can do this better. You can simplify this by introducing extra sets:

`alias (l,i), (m,j); `

set l2(l); l2(l)$(ord(l)<=p) = yes;

set m2(m); m2(m)$(ord(m)<=p) = yes;

set dig2(dig); dig2(dig)$(ord(dig)<=digmax) = yes;

set mapl(l,i); mapl(l2(l),i)$(ord(i)=(1+p*(ord(l)-1))) = yes;

set mapm(m,j); mapm(m2(m),j)$(ord(j)=(1+p*(ord(m)-1))) = yes;

display mapl,mapm;

...

equation obj;

obj.. PROFIT =e= sum((mapl(l,i),mapm(m,j),dig2(dig)) , x(i,j,dig)*ord(dig));

This way you can display the intermediate sets and debug the model more easily. The equation becomes much simpler.

## Tuesday, September 2, 2008

### Privacy.pdf not found

*> You told me to look at*

*http://www.gams.com/docs/privacy.pdf*

*but I get*

*> 403 Forbidden. Please help!!!!*

Sorry. GAMS moved their web pages and this caused problems here and there. I'll email a copy.

## Monday, September 1, 2008

### GAMS bug in redeclaration

`set i;`

set i /1,2/;

set i(*);

The GAMS compiler gets confused if

*i*is used as domain:

` 1 set i;`

2 parameter p(i);

3 set i /1,2/;

4 set i(*);

**** $184

**** 184 Domain list redefined

Work-around: make sure line 1 and 4 use both

*i*or

*i(*)*.