Backend: Integer linear systems are not handled properly

[Martin Sjölund]

The following will generate a (real) linear system of equations. But it contains an integer, which will then truncate the results.

model M
  Real x;
  Integer y;
equation
  2*x + y = 5.0;
  2*x - y = 3.0;
end M;

[Willi Braun]

Dymola says to that model:
{{{The problem is structurally singular for the element type Real.
The number of scalar Real unknown elements are 1.
The number of scalar Real equation elements are 2.
}}}
But how should we treat that issue?
I propose that we just throw a warning that the solution is obtained by truncate some Integer.
Better ideas?

[Martin Sjölund]

As far as I can tell, OMC already links against ​http://lpsolve.sourceforge.net/5.5/ which knows how to solve ILP problems.

I could try to write some C code to see how the solver is used. It seems to support constraints, so we could even use min and max attributes of variables if needed. But we would need a new SimCode structure for these problems (possibly a boolean in SES_LINEAR).

Don't we already support mixed linear systems? Or are those only for boolean variables? We could probably also use lpsolve for these kind of problems.

[Martin Sjölund]

Of course, until we solve add support for it, a warning is sufficient in my opinion

[Willi Braun]

But I mean for that system there is not sufficient solution for y = SOME(Integer). Maybe we shouldn't pass this to the solver like Dymola, because the solver will solve it anyway wrong.

[Martin Sjölund]

Ok, I'll change it to something solvable :)
Does Dymola like this one? We could plug in a solver that can handle this:

/* http://openmodelica.org:8080/cb/issue/1514 model M eq 1+2 */

/* Objective function */
min: ;

/* Constraints */
eq1: +2 X +Y = 5;
eq2: +2 X -Y = 3;

/* Variable bounds */
X >= -Inf;
Y >= -Inf;

/* Integer definitions */
int Y;
Model name: http://openmodelica.org:8080/cb/issue/1514 model M eq 1+2
                 X        Y 
Minimize         0        0 
eq1              2        1  =        5
eq2              2       -1  =        3
Type          Real      Int 
upbo           Inf      Inf 
lowbo         -Inf     -Inf 

Model name:  'http://openmodelica.org:8080/cb/issue/1514 model M eq 1+2' - run #1    
Objective:   Minimize(R0)
 
SUBMITTED
Model size:        2 constraints,       2 variables,            4 non-zeros.
Sets:                                   0 GUB,                  0 SOS.
 
                   1 variables' final bounds.............. RELAXED.
                     OPTIMAL solution found............... 0 
Using DUAL simplex for phase 1 and PRIMAL simplex for phase 2.
The primal and dual simplex pricing strategy set to 'Devex'.
 

Relaxed solution                   0 after          2 iter is B&B base.
 
Feasible solution                  0 after          2 iter,         0 nodes (gap 0.0%)
 
Optimal solution                   0 after          2 iter.

Excellent numeric accuracy ||*|| = 0

 MEMO: lp_solve version 5.5.0.13 for 64 bit OS, with 64 bit REAL variables.
      In the total iteration count 2, 0 (0.0%) were bound flips.
      There were 0 refactorizations, 0 triggered by time and 0 by density.
       ... on average 2.0 major pivots per refactorization.
      The largest [LUSOL v2.2.1.0] fact(B) had 3 NZ entries, 1.0x largest basis.
      The constraint matrix inf-norm is 2, with a dynamic range of 2.
      Time to load data was 0.001 seconds, presolve used 0.000 seconds,
       ... 0.000 seconds in simplex solver, in total 0.001 seconds.
lpsolve OK: 0

Actual values of the variables:
X                               2
Y                               1

[Willi Braun]

No, Dymola even do like this model, with the same answer.
But the point is we can't say without solving (while compile time), if such a system has a solution or not.
So we need to terminate if this system has no suitable solution, not only warn to user.

[Martin Sjölund]

Yes, but even linear systems have this problem (matrix may become singular at any point in time). And nonlinear systems are even worse.

[Willi Braun]

And yes we support such a system.

[Martin Sjölund]

So what's the verdict?
Should we use a normal linear solver like dgesv, truncate and warn the user? Use dgesv, check if the result is integer, else give a runtime error? Use lpsolve, give an error if there is no solution? (2) and (3) are really the same since we have the same number of variables as constraints ((3) is more exact though). Solve only such systems if they have constant jacobian (we can easily solve these during compile-time)?

[Willi Braun]

Yes, you are totally right we need to utilize lp_solve instead of current solution for all kinds of mixed systems.

[Martin Sjölund]

lp_solve is now available as X=System.lpsolve55(A,B,intIndices) (basically solves A*X = B and checks that the correct variables are integers).