Gmail problems

I use gmail almost exclusively and I am very happy with how it works. Recently there surfaced an interesting bug. Attachments can be downloaded individually or combined as a single zip file. It turns out gmail adds a few bytes to the zipped files. After unzipping I see:

C:\Users\erwin\Downloads>dir "AM RFP IFPRI.docx"/s Volume in drive C is HP Volume Serial Number is 8AA3-C6D2 Directory of C:\Users\erwin\Downloads 12/22/2010  01:03 PM           194,667 AM RFP IFPRI.docx                1 File(s)        194,667 bytes Directory of C:\Users\erwin\Downloads\rfpburundi 12/24/2010  12:06 PM           195,072 AM RFP IFPRI.docx                1 File(s)        195,072 bytes      Total Files Listed:                2 File(s)        389,739 bytes                0 Dir(s)  315,865,411,584 bytes free C:\Users\erwin\Downloads>

First one downloaded as single file, second one unzipped from zip. After a few nasty warning messages about corrupted files Word is actually able to open and view the second file also.

See also: http://www.google.com/support/forum/p/gmail/thread?tid=7d4f236a3538ebbc&hl=en

A scheduling problem

For a scheduling that I could not solve in one swoop I tried to develop a small algorithm. Basically: schedule the m<n jobs that come first, then fix the m jobs and schedule the next m jobs. This seems to work with Cplex. The algorithm looks like:

loop(k,    j(j0) = kj(k,j0);    solve sched minimizing makespan using mip;    report(k,'obj') = sched.objval;    report(k,'modelstat') = sched.modelstat;    report(k,'solvestat') = sched.solvestat;    xstart.fx(j,t) = xstart.l(j,t); );

This algorithm will perform the following:

1-10	11-20	21-30	31-40	41-50
solved
fixed	solved
fixed	fixed	solved
fixed	fixed	fixed	solved
fixed	fixed	fixed	fixed	solved

For Cplex we get:

----    180 PARAMETER report             obj   modelstat   solvestat k1    1480.000       1.000       1.000 k2    2940.000       1.000       1.000 k3    4400.000       1.000       1.000 k4    5860.000       1.000       1.000 k5    7320.000       1.000       1.000

However for Gurobi we get:

----    180 PARAMETER report             obj   modelstat   solvestat k1    1480.000       1.000       1.000 k2    2940.000       1.000       1.000 k3    4400.000       1.000       1.000 k4                  19.000       1.000 k5                  19.000       1.000

It turns out that the fixing of a solution causes the next iteration to be infeasible. We can actually isolate step k3 and do:

solve m minimizing makespan using mip;    xstart.fx(j,t) = xstart.l(j,t);    solve m minimizing makespan using mip;

This will make the second model infeasible: i.e. Gurobi does not like its own solution! In this case it is a tolerance question: the final solution is slightly infeasible. We can fix this by:

loop(k,    j(j0) = kj(k,j0);    solve m minimizing makespan using mip;    report(k,'obj') = m.objval;    report(k,'modelstat') = m.modelstat;    report(k,'solvestat') = m.solvestat;    xstart.fx(j,t) = round(xstart.l(j,t)); );

because we know all job step times are in whole seconds. Now both Cplex and Gurobi can solve the problem just fine.

A reasonable good solution looks like:

By accident our algorithmic batch size is the same as the batch size for a good solution (this is not required for the algorithm to work).

The solution is somewhat sensitive where we put the breaks when decomposing the model. This may indicate a rolling horizon algorithm may work better (see http://yetanothermathprogrammingconsultant.blogspot.com/2008/06/rolling-horizon-implementation-in-gams.html).

FTP client

Not many people seem to know this: the Windows file explorer can be used as a simple FTP client by entering a URL of the following format in the  text box where normally the directory is shown:

ftp://login:password@ftp.xxxx.com

Just drag files and whole directories to here and they will be transferred.

Note: when clicking on a URL like that the Internet explorer will start, allowing a “read-only” view of an FTP site.

SQL Injection

From http://www.fsf.org/blogs/sysadmin/savannah-and-www.gnu.org-downtime:

Wed Nov 24 12:59 UTC -- On the evening before Thanksgiving, an IP located in Tbilisi, Georgia started an attack targeting the savannah.nongnu.org website. The perpetrators used SQL injection attacks to download the entire database of usernames and hashed passwords, and we should assume anything else in the Savannah MySQL database.

One way to prevent this is to solely use prepared SQL statements when accessing the DB. This is the strategy I have used in http://yetanothermathprogrammingconsultant.blogspot.com/2010/10/running-large-models-using-web.html. (Other approaches often mentioned are using an appropriate quoting mechanism, e.g. through mysql_real_escape_string and/or disabling multi statements) I browsed through a few recent books on PHP for this project, and I was amazed how much text is devoted to security; it looked like a third of the material seems to be about security. It truly sends shivers down your spine. Programming web applications becomes more and more a question of being paranoid all the time; this is not really a good development as it takes the fun out of things and also reduces productivity (extra work for code that does not implement any new and possibly exciting features for regular, non-hostile users). From a different perspective one could also say that the development tools are still too primitive such that developers are not isolated from these security issues and can concentrate on the real task at hand.

Collaboration tools

For different projects I have been using some of the following document sharing tools:

• DropBox. This allows you to have a directory that is shared between different users. Advantage: on Windows this just looks like a network drive and you can easily copy files from and to the drop box directory. Disadvantages: just files. Also it looks for windows as a local disk so by default dragging means “move” instead of “copy”. Site: http://www.dropbox.com/.
• Wiki-type environments. I have experience with PmWiki (http://www.pmwiki.org) and MediaWiki (http://www.mediawiki.org). Advantage: good support for HTML presentation of rich text. Disadvantage: somewhat complicated to setup and to secure; file upload facilities may be limited depending on the hosting environment.
• Live.com. Microsofts sharing environment. Nice online editing facilities of Word documents.
• My current favorite: sites.google.com. Clean and well-thought through. Has most of the facilities I need. Easily shared with just the people you want to share it with. Limits: 20 GB files, and 100 GB per site; I had to park some of the larger data files somewhere else.

Reorganizing a large parameter in GAMS

To export some fairly substantial data in a certain format I wanted to form:

DownLoadData(downloadfiles,SP_Varname,SP_SubVarName,SP_SubSubVarName,SP_Indicator,                                 SP_Description,SP_Source,SP_Unit,cty,Fwdy) =       sum((mapdownloadfiles(downloadfiles,SP_control),SP_Map,SP_Graph,SP_Options,SP_Type,mapy(FwdY,SP_year)),           cube(SP_Control,                SP_VarName,                SP_SubVarName,                SP_SubSubVarName,                SP_Year,                SP_Indicator,                SP_Source,                SP_Description,                SP_Unit,                SP_Map,                SP_Graph,                SP_Options,                SP_Type,                cty                ) );

The original Cube parameter has 326,755 elements, and this operation basically reshuffles things into a lower-dimensional parameter with some mappings applied. Upon completion the new parameter Download_Data has the same number of elements: 326,755. Unfortunately this operation was very expensive, so much so that I never saw it finish.

The following code does exactly the same in just 0.5 seconds:

loop((mapdownloadfiles(downloadfiles,SP_control),mapy(Download_Year,SP_year)),    Download_Data(downloadfiles,SP_Varname,SP_SubVarName,SP_SubSubVarName,cty,SP_Indicator,                              SP_Description,SP_Source,SP_Unit,Download_Year) =         sum((SP_Map,SP_Graph,SP_Options,SP_Type),           cube(SP_Control,                SP_VarName,                SP_SubVarName,                SP_SubSubVarName,                SP_Year,                SP_Indicator,                SP_Source,                SP_Description,                SP_Unit,                SP_Map,                SP_Graph,                SP_Options,                SP_Type,                cty                ) ); );

GAMS has some built-in facilities to automatically rearrange assignments like these to achieve better performance, but in this case it was not very effective. So human intervention was needed to optimize this assignment.

Amazon not just selling books

Dear Amazon EC2 Customer, We are excited to announce the immediate availability of Cluster GPU Instances for Amazon EC2, a new instance type designed to deliver the power of GPU processing in the cloud. GPUs are increasingly being used to accelerate the performance of many general purpose computing problems. However, for many organizations, GPU processing has been out of reach due to the unique infrastructural challenges and high cost of the technology. Amazon Cluster GPU Instances remove this barrier by providing developers and businesses immediate access to the highly tuned compute performance of GPUs with no upfront investment or long-term commitment. Amazon Cluster GPU Instances provide 22 GB of memory, 33.5 EC2 Compute Units, and utilize the Amazon EC2 Cluster network, which provides high throughput and low latency for High Performance Computing (HPC) and data intensive applications. Each GPU instance features two NVIDIA Tesla® M2050 GPUs, delivering peak performance of more than one trillion double-precision FLOPS. Many workloads can be greatly accelerated by taking advantage of the parallel processing power of hundreds of cores in the new GPU instances. Many industries including oil and gas exploration, graphics rendering and engineering design are using GPU processors to improve the performance of their critical applications. Amazon Cluster GPU Instances extend the options for running HPC workloads in the AWS cloud. Cluster Compute Instances, launched earlier this year, provide the ability to create clusters of instances connected by a low latency, high throughput network. Cluster GPU Instances give customers with HPC workloads an additional option to further customize their high performance clusters in the cloud. For those customers who have applications that can benefit from the parallel computing power of GPUs, Amazon Cluster GPU Instances can often lead to even further efficiency gains over what can be achieved with traditional processors. By leveraging both instance types, HPC customers can tailor their compute cluster to best meet the performance needs of their workloads. For more information on HPC capabilities provided by Amazon EC2, visit aws.amazon.com/ec2/hpc-applications. Learn more about the new Cluster GPU instances for Amazon EC2 and their use in running HPC applications. Sincerely, The Amazon EC2 Team

Unfortunately I don’t have any direct use for this (it would be fun to try this out on some problems). It is interesting that Amazon is so active in the HPC arena.

U shaped penalties

For a large LP model a customer of mine wanted to have instead of a V shaped penalty around a goal, something that approximates a U shaped penalty. With standard slacks:

 

we essentially have a V shaped penalty:

We can approximate a more U shaped form by introducing extra slacks with a smaller penalty:

with unit penalty q smaller than p, e.q. q=0.1*p. This will give:

The idea is that the algorithm will select the cheap slacks t first, before spending on the expensive slacks s. This will keep the LP linear and we don’t need extra binary variables. For q=0 we essentially allow small deviations for free, but start paying when we go beyond the threshold T.

This approach can easily be extended to more piecewise linear segments (as long as we keep things convex). Also the range T for the slacks t may be generalized to:

StatPlanet

When demoing http://www.sacmeq.org/statplanet/StatPlanet.html (for which I developed a GAMS link), I reminded to Hans Rosling. Check here for some fantastic presentations: http://www.gapminder.org/videos/. The software has been acquired by Google. Here is an example:

Amazon ec2

Amazon “Elastic Compute Cloud” has some big-ass machines available: http://aws.amazon.com/ec2/#instance. See also the rates on that page.

Running large models using a web interface

For a project I am considering a fairly large GAMS models where we want end-users to run certain scenarios. Of course there must be HTML GUI to allow the user to specify the input for each scenario to run. That is fairly standard stuff – once you know what type of parameters you want to expose to the user to change. The bigger issue is that the model is too large to run immediately if the user presses the Submit button. Rather we need to create a small queueing mechanism where jobs are waiting to run. After the model has been solved the user can inspect the results.
Running background jobs from a web-server environment is a little bit more exotic, although other examples exist (see e.g. NEOS and this effort). For this project we use a Windows IIS web server with PHP.

Architecture

A background solve process can look like:

1. If background processes are already running, stop
2. Start a background process as follows:
3. While there are waiting jobs
1. select a job where status=waiting
2. set status=running
3. solve model
4. store results
5. set status=solved

We run this thing after every INSERT in the job table.
Actually we can have several of those processes if the server is beefy enough to allow for parallel solution of different models. This does not only depends on the number of cores available but also memory considerations play a role. Essentially line 1 becomes:

1. If number of background processes ≥ maxnumproc, stop

Speed of loops in GAMS

In general loops are slow in GAMS. See http://yetanothermathprogrammingconsultant.blogspot.com/2009/08/gams-speed-of-loops.html. There are some exceptions however. Here is an example:

Note that map(i,k) has a diagonal structure.

The timings are:

method 1: loop	0.156 seconds
method 2: sum	3.495 seconds

Advanced basis/starting point in GAMS

Savepoint/loadpoint is an easy way to use an initial point for an NLP or advanced basis for an LP. The commands are:

 

The first time you run this the solution is saved in the GDX file. The next time around, this GDX is used to set as many levels and marginals that can be matched by name. For this very large MCP model it turned out to be useful to bring down solution times of related experiments.

Sorry for the ugly syntax. I would have preferred something more symmetric for reading and writing.

Assertive student

Hi Erwin,

In short, I'm student and wish not to receive "ask your supervisor" answer.

Is there a function, like the error function, to compute the probability of random variable that has a beta distribution? if not, is there any way to statistically represent a bounded variable?

Thanks for you consideration,

This is probably related to http://yetanothermathprogrammingconsultant.blogspot.com/2010/09/student-mails.html. My hint in this case was: See http://www.amsterdamoptimization.com/pdf/specfun.pdf.

Dissemination of model results

After spending time and effort on a model, it may well be worth to pay attention to presentation and dissemination of the model results. Currently I am working on a project where we can export results for a very large and complex GAMS model to an Adobe Flash application called StatPlanet. Obviously a web-based tool to inspect results has lots of advantages (nothing to install, available from everywhere).

GDX2StatPlanet

GAMS/GDX

GAMS has a binary data file format for which an API is available. In the last few weeks I encountered a number of similar issues with this in different projects.
The interface does require more knowledge of the internals of GAMS than most people have. An example is here. When we want to read a GAMS parameter A(I,J) into a C/Fortran/Delphi matrix A[1..m,1..n] we see this is not completely trivial. The basic obstacle is that GAMS has strings as indices, while most programming languages have simple integers as indices. I.e. we need to map J = /new-york, chicago, topeka/ into j=1..3. For large problems that means some kind of hashing or other way to quickly look up strings. There is also a “raw” mode interface that does not give strings but rather internal GAMS numbers. So we need to map J = /1001,1002,1005/ to j=1..3. This again is not a completely trivial task. If you are familiar with the internals of GAMS one can know that GAMS internal numbers (so-called UEL numbers) are by definition increasing. So one could use a binary search to quickly find numbers. A fortran routine for this can look like:

Another possibility would by to have an array of the size of MAXUEL that has mapping information from UEL-> j=1..3. That requires more space but has quick look up. If the GDX file contains only one set, one could use socalled “mapped” mode, but that cannot be used in general.

It is also worth noting that GAMS only stores and exports nonzero elements. So we cannot assume all A(i,j)’s are in the GDX file.

So to recap, my suggested approach to read a matrix A(i,j) from a GDX file is:

1. Read 1D set I (in raw mode) and store UEL’s in an integer array IA
2. Read 1D set J (in raw mode) and store UEL’s in an integer array JA
3. Allocate array A(i,j) and set all elements to zero
4. Read 2D parameter A (in raw mode). For each record do:
   1. Retrieve (ueli, uelj, aij)
   2. i := bsearch(IA, length(IA), ueli)
   3. j := bsearch(JA, length(JA), uelj)
   4. A[i,j] := aij

If you don’t want to use binary search an alternative algorithm would be:

1. Retrieve MAXUEL (gdxSystemInfo)
2. Allocate integer arrays IA and JA of length MAXUEL
3. Read 1D set I (in raw mode). For each record do:
   1. Retrieve uel
   2. IA[uel] := recordnumber
4. Read 1D set J into array JA as described in step 3
5. Read 2D parameter A (in raw mode). For each record do:
   1. Retrieve (ueli, uelj, aij)
   2. i := IA[ueli]
   3. j := JA[uelj]
   4. A[i,j] := aij

Notes:

• The generated interface files are quite ugly, have strange design and implementation and have no comments on usage. The code in the provided files just looks very awkward. Below are a few more specific issues illustrating this.
• The filenames are incomprehensible: gdxdocpdef.pas, gdxdcpdef.pas ??? 
• The functions should return boolean/logical values instead of integers.
• The previous point is the more crucial as some routines have 0=success and others have a return code <>0 meaning success. That throws me off balance quite often.
• Fortran support for many compilers is missing.
• The examples are too simplistic. The above is not explained at all.
• The generated interface contains some code that I would refuse from my students when I was teaching programming at undergraduate level. An example is:
     n: ShortString
  A symbol n should always be an integer, never a string (or a floating point number).
• The Delphi interface does not compile under the latest compiler.

Delphi 2010

GAMS users routinely generate very large GDX files with millions of records. To get data in and out of these files efficiently still requires substantial skills and effort.

Large NLP/MCP models

For a very large large MCP model (a system of nonlinear equations with complementarity conditions) we observed it became unsolvable, even with a very good starting point. The reason was that the model was kept as small as possible in terms of equations and variables. As a result the model was much more dense and nonlinear (in terms of nonlinear nonzero elements) than needed. Basically, modeling systems and solvers that exploit sparsity should be fed large, sparse models instead of corresponding smaller, denser models. After substituting out some common expressions into separate equations (hence increasing the size in terms of equations and variables) we could solve the larger models (≈400K variables and equations) fairly easily.

An linear example of such a problem is shown here: http://yetanothermathprogrammingconsultant.blogspot.com/2010/04/l1-portfolio-formulation.html.

An example of a nonlinear construct that should raise a red flag is something like:

In many cases it is better to write this as:

This will keep the number of nonlinear nonzero elements smaller at the expense of extra variables and equations.

In the model I was looking at a combination of these two situations was implemented: a large linear expression inside a nonlinear expression, and this nonlinear expression was repeated several times in the model. This caused both extra nonzero elements, and they appeared also in very nonlinear fashion.

For a smaller version of the model, the differences in statistics is striking (I was unable to get the large version to solve with the original formulation):

old model

MODEL STATISTICS BLOCKS OF EQUATIONS          25     SINGLE EQUATIONS      106,219 BLOCKS OF VARIABLES          25     SINGLE VARIABLES      106,236 NON ZERO ELEMENTS     5,238,688     NON LINEAR N-Z      5,034,974 DERIVATIVE POOL           1,828     CONSTANT POOL         104,745 CODE LENGTH         120,759,139 GENERATION TIME      =       52.151 SECONDS EXECUTION TIME       =       52.463 SECONDS RESOURCE USAGE, LIMIT        281.806      1000.000

new formulation

MODEL STATISTICS BLOCKS OF EQUATIONS          26     SINGLE EQUATIONS      107,605 BLOCKS OF VARIABLES          26     SINGLE VARIABLES      107,622 NON ZERO ELEMENTS       768,295     NON LINEAR N-Z        561,809 DERIVATIVE POOL           1,859     CONSTANT POOL         104,745 CODE LENGTH          14,348,068 GENERATION TIME      =        7.566 SECONDS EXECUTION TIME       =        7.972 SECONDS RESOURCE USAGE, LIMIT         57.216      1000.000

Unfortunately solvers and modeling systems don’t give much feedback on these issues. Such a diagnoses requires some modeling knowledge and experience (one reason why I am sometimes hired to help out).

Student mails

I get many mails from students asking for help. I usually try to give a small hint and direct them to get more help from their supervisor. Here are some examples.

I work for a supply chain management company, and we need some help on the attached exercise we have been given as an internal training exercise. Can you help? Attached is the project and all of the data and mod files. The deadline is Friday. If you can help, please provide a price quotation.

Answer: I would suggest to takes this up with your instructor. I believe it can be conceived as inappropriate if I would step into this.

After my standard reply of "Talk to your supervisor. He/she is actually getting paid to help you." one sometimes get long stories of poor supervision by professors:

I am really, really sorry to bother you, but my supervisor does not know how to use CPlex, only one lecturer in Math department knows how to use CPlex, but he is very busy at the moment (his wife just got cancer), you are the only expert that I can count on, it is a critical part of my PhD thesis, I really, really need your help for my problem, would you please have a look at my question and give me some idea? Thank you billion times! Your help will be greatly appreciated!

Or sometimes shorter:

Me and my wife are in very difficult situation. She runs out of the scholarship and I have to go back to Kuwait with her and the family.

Of course if I send back a file with some results, I really should hide all details about who actually did the work. Secrecy is important here: 

Also can you change the following to STUDENTSNAME Licensee: Erwin Kalvelagen                               G080731/0001CJ-WIN Nobody know you are helping me please. My future is between your hands Can you run GAMS economical model / Spread sheet please and put my name STUDENTSNAME everywhere please?

Some people are quite brazen in their request. This guy wants me to do his complete PhD work and he will find someone else to help me if I am too expensive:

What is your lump thump price tell I get my PhD to get two publication accepted and modify my GAMS RESULTS tell the committee is satisfied and teaching me what you did and prepare me also for defence? What is the budget for the first publication is $? while how much you should be able to write the second publication for $?. you will make sure that those publications are of such a high standard that you will be able to find a journal which will accept both papers. In case that this is more than my budget allows, then we can just do the GAMS part for me. Once this is done we might as well find someone else who can write you the two publications.

Does not look PhD supervisors are doing a good job if students think they can get away with such requests!

GAMS macro problems

I am working on a model that is using the preprocessor tool described in http://www.mpsge.org/inclib/gams-f.htm. There were a few minor issues related to this tool:

1. It does not really work very well with the GAMS IDE. You cannot click on the red error messages to bring you to the input line with the error.
2. It does not handle combined use in equations and parameter assignments in post-solve reporting very well.

The documentation of the new $macro facility in GAMS (http://www.gams.com/docs/release/release.htm, section “The GAMS Macro Facility” ) seems to indicate that it can be used to replace the preprocessor tool: just use $macro and use set.local where needed. Unfortunately that is not completely true. Especially the sentence “The .local feature has been added to ensure local use and eliminates the need for new alias definitions.” is plain wrong.

Here is a simple model to illustrate the problem:

$sysinclude gams-f set i /i1*i3/; alias (i,j); parameter x(i,j); x(i,j) = normal(0,1); *----------------------------------------------------- parameter p(i) 'without macro/preproc'; p(i) = sum(j, x(i,j)); display p; *----------------------------------------------------- y(i) == sum(j, x(i,j)); parameter p1(j) 'with preproc'; p1(j) = y(j); display p1; *----------------------------------------------------- $macro y2(i) [sum(j.local, x(i,j))] parameter p2(j) 'with macro'; p2(j) = y2(j); display p2;

This shows:

----     26 PARAMETER p  without macro/preproc i1  0.478,    i2 -3.533,    i3 -0.219 ----     49 PARAMETER p1  with preproc i1  0.478,    i2 -3.533,    i3 -0.219 ----     57 PARAMETER p2  with macro i1 -1.804,    i2 -1.804,    i3 -1.804

The cause can be found in the listing file:

35  ALIAS(FCN000,j); 36   37   38   39  parameter p1(j) 'with preproc'; 40  p1(j) = 41   42  *       #### Preprocessor function replacement for Y #### 43   44  ( 45   sum(FCN000, x(j,FCN000)) 46  ) ……. 55  parameter p2(j) 'with macro'; 56  p2(j) = [sum(j.local, x(j,j))];

Clearly the GAMS-F preprocessor handles this case correctly by introducing an alias. The $macro/.local combination does not do this and gives the wrong results. Instead of this .local hack it is needed to manually introduce an alias and use the aliased set where needed.

It is noted that overuse of these tools can lead to models that appear smaller (fewer equations and variables). However those models may be much denser and much more nonlinear (say in terms of nonlinear nonzero elements). Often larger, sparser and less non-linear models solve easier.

Google or-tools

The first OR-tool from Google is a preliminary constraint programming solver, available from here: https://code.google.com/p/or-tools/. The name “or-tools” seems to indicate we can expect more to come.

Efficient frontier-like analysis

I needed to perform an optimization run of the form:

for different λ’s. Essentially something like an efficient frontier but in this case we have a MINLP (with semi-integer variables), so no smooth curve but instead some points.

In this case f(x) and g(x) were having different units and different magnitudes. In that case it is better to solve first:

and then use:

This will get a (hopefully) better spread of the λ’s. See: http://yetanothermathprogrammingconsultant.blogspot.com/2010/02/efficient-frontiers-in-mip-models.html.

Of course if λ=0 or λ=1 is part of the λ’s we need to trace, then this actually does not require extra models to solve, just a different ordering: first do λ=0 and λ=1 to get f₀ and g₀ and then do the λ’s in between. As the λ’s are user specified, we needed a little bit of logic to handle this.

Long running horizon

I am experimenting with a timber model. These models have typically long horizons, which can cause some problems with some solvers. As the model has a nonlinear objective I tried a few NLP solvers:

Solver	Objective (maximize)
Conopt	104599.7740
Minos	53746.0234
Snopt	6651.1413

This looks pretty bad. Fortunately, a simple option setting: SCALE OPTION = 0, can make things much better:

Solver	Objective (maximize)
Minos (no scaling)	104599.7740
Snopt (no scaling)	104599.6879

Intuitively this makes some sense: MINOS and SNOPT scale by default only the linear variables. This strategy is actually not without merit: the Jacobian elements for these nonlinear variables can change quite a lot during the optimization, so scaling based on the initial Jacobian elements may not be very wise. However, in some cases making this distinction between linear and nonlinear variables can make things worse. In this case it is better not to scale at all.

I was first tinkering with the optimality tolerance, as I have seen with other forestry models (LP and NLP) that for these ridiculously long time horizons, tightening the optimality tolerance can help getting better solutions. For this model that was not the case.

Cplex new version

In some cases a new version of a solver can perform worse on some models. Here is a particular unfortunate example, where Cplex 12.1 finds several very good solutions in the initial heuristics, but Cplex 12.2 just fails to find any integer solutions within the allotted time: 

ILOG CPLEX       Nov  1, 2009 23.3.3 WEX 13908.15043 WEI x86_64/MS Windows Cplex 12.1.0, GAMS Link 34 Cplex licensed for 1 use of parallel lp, qp, mip and barrier. Cplex MIP uses 1 of 4 parallel threads. Change default with option THREADS. Reading data... Starting Cplex.… … Root relaxation solution time =  102.77 sec.         Nodes                                         Cuts/    Node  Left     Objective  IInf  Best Integer     Best Node    ItCnt     Gap       0     0    15053.5468   150                  15053.5468   151633         *     0+    0                       150008.6529    15053.5468   151633   89.96%       0     0    15065.0166   146   150008.6529     Cuts: 479   158521   89.96% *     0+    0                        15387.9888    15065.0166   158521    2.10%       0     0    15065.1974   152    15387.9888 Flowcuts: 146   159618    2.10%       0     0    15065.2096   146    15387.9888      Cuts: 70   159894    2.10%       0     0    15065.2104   146    15387.9888   Flowcuts: 9   159952    2.10% *     0+    0                        15175.4264    15065.2104   159952    0.73% Flow cuts applied:  472 Mixed integer rounding cuts applied:  4 Gomory fractional cuts applied:  7 MIP status(102): integer optimal, tolerance

 

IBM ILOG CPLEX   Jul  4, 2010 23.5.1 WEX 18414.18495 WEI x86_64/MS Windows Cplex 12.2.0.0, GAMS Link 34 GAMS/Cplex licensed for continuous and discrete problems. Cplex MIP uses 1 of 4 parallel threads. Change default with option THREADS. Reading data... Starting Cplex... … Root relaxation solution time =  107.14 sec.         Nodes                                         Cuts/    Node  Left     Objective  IInf  Best Integer     Best Node    ItCnt     Gap       0     0    15053.5468   111                  15053.5468   120941               0     0    15066.3609   105                   Cuts: 787   126925               0     0    15067.4221    99                     Cuts: 7   127524               0     0    15067.4512   100                     Cuts: 8   127557               0     0    15067.4512   100                  MIRcuts: 2   127559         Heuristic still looking. Heuristic still looking.       0     2    15067.4512   100                  15067.4561   127559         Elapsed real time = 237.15 sec. (tree size =  0.00 MB, solutions = 0)       1     3    15072.3309    97                  15067.4561   131486               2     4    15068.7687    75                  15068.8163   134032               3     5    15069.6752    84                  15068.8163   134657               4     6    15071.3502    88                  15069.6828   135922               5     7    15072.3709    75                  15069.6828   136786               6     8    15086.2295    77                  15069.6828   145563               7     9    15072.9117    71                  15069.6828   146885               8    10    15073.9185    79                  15069.6828   148355               9    11    15081.7555    70                  15069.6828   154282              12    14    15074.6296    72                  15069.6828   154696         Elapsed real time = 314.54 sec. (tree size =  2.27 MB, solutions = 0)      13    15    15088.4559    66                  15069.6828   159070         … MIP status(108): time limit exceeded, no integer solution CPLEX Error  1217: No solution exists. Resource limit exceeded, no integer solution found.

Such behavior is always difficult to explain to a customer. The argument that on average a new version tends to do better is of little consolation.

PS. The time limit was very tight on this model: 400 seconds (extremely tight considering the LP takes 100 seconds). This was no problem with previous versions of Cplex as the initial heuristic always found solutions very quickly. These solutions were so good they already obeyed the requested gap of 0.8%. Cplex 12.2 finds the solution relatively quickly in the B&B after about 500 seconds on my machine with 1 thread (the client on his machine obtained no solution in 600 seconds with 2 threads and some other solver options; I assume a little bit more time may help here).

PS2. All my runs were with default settings (except 1 thread, gap=0.8%, time limit=400 seconds).

QP Models with GAMS/Knitro

> My large QP portfolio model does not solve with GAMS/KNITRO.

Yes, for large problems you will get:

Memory estimate for Hessian is too small.
Increase the <modelname>.workfactor option from the current value of 1

For 2000 stocks with WORKFACTOR=10 I got the same error.

Even if the model is formulated as a QCP model, the GAMS link is not smart enough in setting up the hessian efficiently. Instead it looks like the QP is handled as a general NLP.

There are two easy workarounds. First of all you can solve this with IPOPT. The GAMS link for IPOPT is smarter in dealing with second derivatives than the link for KNITRO (as of GAMS 23.5). Note that it is better to solve as a QCP model than as an NLP model: GAMS is not smart enough to recognize this as a QP and exploit that:

model type	GAMS generation time	IPOPT time	function evaluation time	Total solver time
NLP	38	21.8	39.1	61.0
QCP	38	21.6	9.8	31.4

A second possibility is to reformulate model (1):

into model (2):

This formulation can be solved with KNITRO without the problems with the hessian storage.

---

As an aside we can improve formulation 1 by replacing the variance-covariance matrix Q by

We exploit here symmetry in the quadratic form x’Qx. This exercise substantially reduces the GAMS generation time, and also further shaves off some time from the solver:

model type	GAMS generation time	IPOPT time	function evaluation time	Total solver time
QCP –triangular Q	18.8	21.3	4.6	26

We can also apply the triangular storage scheme to model (2). With a solver like SNOPT this leads to a very fast total solution time (i.e. generation time + solver time).

model type	GAMS generation time	SNOPT time
NLP –triangular Q in linear constraints	1	5

If needed we can even speed this up by not using the covariances but rather the mean adjusted returns directly:

I used T=300 time periods. This formulation yields:

model type	GAMS generation time	SNOPT time
NLP – Mean Adj Returns	0.5	2.3

Even for the simple QP portfolio model, careful modeling can help to improve the performance in a significant way.