linear programming

Librairie C++ de calcul formel/ C++ symbolic computation library

Modérateur : xcasadmin

lukamar
Messages : 331
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

linear programming

Message par lukamar » mar. août 08, 2017 8:35 pm

Hi Bernard,

I would like to contribute to the Giac development by providing the c++ code for function 'lpsolve' that emulates (almost fully) the function LPSolve from Maple's Optimization package. From the comment in misc.cc source file I see that an implementation of the two-phase simplex method is on the TODO list. I implemented it (upon already written 'simplex_reduce') alongside with the branch&bound method for solving (mixed) integer and binary LP problems. These are the main tools for 'lpsolve' (no active set or interior point methods were implemented).

With this addition (few hundreds lines of code), Giac would be able to solve great variety of practical problems, those which can be modeled as (integer) LP problems, Also, the compatibility with Maple would improve. As I am particulary interested in Giac having an usable LP solver, I decided to write one myself. I would be happy if you'd find time to review the code (it's GNU/GPL) and tell me if it qualifies for the testing version of Giac. Also, how do I submit the code? Attaching documents to messages in this forum seems to be disabled.

parisse
Messages : 5734
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

Re: linear programming

Message par parisse » mer. août 09, 2017 6:10 am

That would be really great!
I think you can attach zip files. Or enter the source in a message inside code delimiters, or email it to me!

lukamar
Messages : 331
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

Re: linear programming

Message par lukamar » lun. août 14, 2017 5:16 pm

Here is the code. Sorry for the delay, I had to fix some bugs.
The enclosed Makefile creates executable a.out which can be run to test the function lpsolve. Just enter the problem as "(obj,constr,....)". There are several examples of LP problems included, mostly from the Maplesoft web page explaining their function LPSolve. They should all work except the last, eleventh example in the documentation entry for function lpsolve, I have no idea why it does not work.
There are couple of things I should mention regarding lpsolve.cc. As the inclusion of 'plus_inf' and 'minus_inf' constants in the code made it segfault for some reason, I used a workaround by adding function 'make_infinity' which creates +infinity symbol. It should be removed and all instances of 'inifnity' outside the comments replaced by 'plus_inf', as in other Giac functions. Also, there was a bug in function 'simplex_reduce' which I had to correct, hence the fixed version 'simplex_reduce2'. When determining basis after pivoting, the 'simplex_reduce' only checked for zeros in an idn column, but not that the nonzero element is equal to one. I added that check (see line 171 in lpsolve.cc), which should be included in the code for 'simplex_reduce'.
For example, to solve the LP problem included in cascmd_fr (when explaining 'simplex_reduce' function), one should run a.out and enter:

Code : Tout sélectionner

(2x+y-z+4,[x<=1,y>=2,x+3y-z=2,2x-y+z<=8,-x+y<=5])
Also, the function 'lpsolve' requires a few new option keywords to be added to Giac. These are listed in lpsolve.h. I used the available option 'integer', but other options, such as 'nonnegative' and so on, should be defined in similar fashion. Since they're missing, use the corresponding integer values defined in lpsolve.h. For example, to specify 'assume=nonnegative' one should type '0=7'.
Pièces jointes
lpsolve.zip
(15.56 Kio) Téléchargé 229 fois

parisse
Messages : 5734
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

Re: linear programming

Message par parisse » mar. août 15, 2017 1:58 pm

I have mirrored your bugfix in misc.cc, so that you can remove simplex_reduce2.
For the lp keywords, I can add them to the lexer, except that:
* assume is already a giac keyword, you can't parse it as an integer, but it's not difficult to check for it.
* nonnegint. already exists, value is _NONNEGINT=4*256+2,
I can add the lp specific integer constants to maple_enum in dispatch.h, values may therefore change with respect to lpsolve.h (lpsolve.cc unchanged).
This could bring something like:

Code : Tout sélectionner

      if (R==at_equal) {
        // parse the argument in form "option=value"
        vecteur ops(*((*it)._SYMBptr->feuille._VECTptr));
        if (ops[0].type==_IDNT && intrv2vec(ops[1],lr,contextptr)) {
          // the boundaries for variable ops[0] are set
          if (!is_zero(lr[0]) && is_strictly_greater(lr[0],-infinity,contextptr))
            bconstr.push_back(symbolic(at_superieur_egal,makevecteur(ops[0],lr[0])));
          if (is_strictly_greater(infinity,lr[1],contextptr))
            bconstr.push_back(symbolic(at_inferieur_egal,makevecteur(ops[0],lr[1])));
          if (is_positive(lr[0],contextptr))
            nvars.push_back(ops[0]);
          continue;
        }
        if (ops[0]==at_assume){ // put here the case LP_ASSUME:
          ...
         continue;
        }
        if (is_integer(ops[0])){
          ...
        }
     }
Once all is fixed, I can simply add lpsolve.cc/lpsolve.h to Makefile.am and lpsolve should be available in Giac/Xcas. I don't know why you get segfault with the builtin inf, perhaps a loader order bug, anyway this should disappear once the code is integrated in the library.

parisse
Messages : 5734
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

Re: linear programming

Message par parisse » mar. août 15, 2017 2:21 pm

Actually, I think it's safer to add the prefix lp_ to all constants that the user will key in (while accepting both lp_integer and integer or lp_nonnegint and nonnegint should be supported). Code added in lexer_tab_int.h:

Code : Tout sélectionner

      {"lp_binary"               ,0, _LP_BINARY, _INT_MAPLECONVERSION, T_TYPE_ID},
      {"lp_binaryvariables"               ,0, _LP_BINARYVARIABLES, _INT_MAPLECONVERSION, T_TYPE_ID},
      {"lp_depthlimit"               ,0, _LP_DEPTHLIMIT, _INT_MAPLECONVERSION, T_TYPE_ID},
      {"lp_integer"               ,0, _LP_INTEGER, _INT_MAPLECONVERSION, T_TYPE_ID},
      {"lp_integervariables"               ,0, _LP_INTEGERVARIABLES, _INT_MAPLECONVERSION, T_TYPE_ID},
      {"lp_maximize"               ,0, _LP_MAXIMIZE, _INT_MAPLECONVERSION, T_TYPE_ID},
      {"lp_nonnegative"               ,0, _LP_NONNEGATIVE, _INT_MAPLECONVERSION, T_TYPE_ID},
      {"lp_nonnegint"               ,0, _LP_NONNEGINT, _INT_MAPLECONVERSION, T_TYPE_ID},
      {"lp_variables"               ,0, _LP_VARIABLES, _INT_MAPLECONVERSION, T_TYPE_ID},
In dispatch.h:

Code : Tout sélectionner

  enum maple_conversion {
    _TRIG=100,
    _EXPLN=101,
    _PARFRAC=102,
    _FULLPARFRAC=103,
    _BASE=104,
    _CONFRAC=105,
    _MAPLE_LIST=256,
    _POSINT=1*256+2,
    _NEGINT=2*256+2,
    _NONPOSINT=3*256+2,
    _NONNEGINT=4*256+2,
    _LP_BINARY=106,           //   * binary
    _LP_BINARYVARIABLES=107,  //   * binaryvariables
    _LP_DEPTHLIMIT=108,       //   * depthlimit
    _LP_INTEGER=109,          //   * integer
    _LP_INTEGERVARIABLES=110, //   * integervariables
    _LP_MAXIMIZE=111,         //   * maximize
    _LP_NONNEGATIVE=112,      //   * nonnegative
    _LP_NONNEGINT=113,        //   * nonnegint
    _LP_VARIABLES=114,        //   * variables
  };
Code also added in gen.cc in printint32.

Code : Tout sélectionner

    ...
      case _NONNEGINT:
	return "nonnegint";
      case _LP_BINARY:
	return "lp_binary";
      case _LP_BINARYVARIABLES:
	return "lp_binaryvariables";
      case _LP_DEPTHLIMIT:
	return "lp_depthlimit";
      case _LP_MAXIMIZE:
	return "lp_maximize";
      case _LP_NONNEGATIVE:
	return "lp_nonnegative";
      case _LP_NONNEGINT:
	return "lp_nonnegint";
      case _LP_VARIABLES:
	return "lp_variables";
      case _LP_INTEGER:
	return "lp_integer";
      case _LP_INTEGERVARIABLES:
	return "lp_integervariables";
      }

lukamar
Messages : 331
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

Re: linear programming

Message par lukamar » mer. août 16, 2017 2:32 pm

Here is the final version of the code. I fixed the bug which caused wrong result to the example 11, it was in branch_and_bound function when comparing with the incumbent solution. Now all examples should work. I put plus_inf and minus_inf on the respective places and the code compiles fine. Also, i made a few improvements, did some refactoring, added more examples and updated the documentation blocks.
I also plan to write an entry to cascmd TeX file (english version) about the function 'lpsolve'. I'll send it to you when I finish, most probably in the next day or two.

Edit: now I see that I forgot to change the names in lpsolve.h. lpmax should be simplex_twophase, lpmin should be simplex_dual and ilpbb should be branch_and_bound.
Pièces jointes
lpsolve-fixed.zip
(13.83 Kio) Téléchargé 233 fois

parisse
Messages : 5734
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

Re: linear programming

Message par parisse » dim. août 20, 2017 5:57 pm

It is now included in the javascript version if you want to check:
https://www-fourier.ujf-grenoble.fr/~pa ... casfr.html
I will update giac source and version during the first week of September (perhaps before if I can).

lukamar
Messages : 331
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

Re: linear programming

Message par lukamar » dim. août 20, 2017 7:32 pm

Thanks. But the web interface does not recognize the 'lpsolve' command. I cleared the cached web content and reloaded, but with no effect. How do I update CAS? I'm using Firefox.

parisse
Messages : 5734
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

Re: linear programming

Message par parisse » lun. août 21, 2017 5:51 am

The UI does not display lpsolve as a command (I must add it in a file somewhere) but it should work.

lukamar
Messages : 331
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

Re: linear programming

Message par lukamar » mar. août 22, 2017 11:07 am

Yes, it works indeed.
I have one small syntactical improvement for lpsolve. After the line 1022 it should be inserted

Code : Tout sélectionner

vector<int> vc=findvars(cvars,vars);
for (vector<int>::const_iterator it=vc.begin();it!=vc.end();++it) {
  vartype[*it]=0;
}
Namely, cvars is not being actually used anywhere (it holds a list of continuous variables), which is not a problem since all variables are continuous by default. But with this additional code the option "lp_variables=[...]" may be used to override already specified integer variables, which can be useful when specifying mixed problems. For example, a problem with 100 integer variables and 2 continuous variables may be specified with options "assume=integer" and "lp_variables=[cont_var1,cont_var2]".

Also, I noticed that lpsolve returns +infinity for unbounded cases of minimization problems. The correct value should be -infinity. That has troubled me before but I thought that assigning minus_inf to soln[0] (see line 1048) is OK. How to fix it? (I could not test it on my computer because of the mysterious segfaulting when using plus_inf, minus_inf, zero etc., which I already reported and hadn't been able to overcome since.)

Finally, I think that the message "branch and bound: ... tree node(s) examined" at line 1063 would be better without the word "tree".

lukamar
Messages : 331
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

Re: linear programming

Message par lukamar » mar. août 22, 2017 1:35 pm

I finished the entry about 'lpsolve' for cascmd_en manual. Below is the TeX file containing the updated section "Linear Programmation". I also added the label on line 47 (subsection: simplex_algorithm) to make a reference to an existing example. The new subsection "Solving general linear programming problems" is appended to the section.

Code : Tout sélectionner

\section{Linear Programmation}\index{simplex\_reduce|textbf}
Linear programming problems are maximization problem of a linear
functional under linear equality or inequality constraints.
The most simple case can be solved directly by the so-called simplex
algorithm. Most cases requires to solve an auxiliary linear
programming problem to find an initial vertex for the simplex
algorithm.

\subsection{Simplex algorithm: {\tt simplex\_reduce}}
{\bf The simple case}\\
The function {\tt simplex\_reduce} makes the reduction 
by the simplex algorithm to find : 
\[ \mbox{max}(c.x), \quad  A.x \leq b,\ x \geq 0,\ b\geq 0 \]
where $c,x$ are vectors of $\mathbb R^n$, $b\geq 0$ is a vector of 
$\mathbb R^p$ and $A$ is a matrix of $p$ rows and $n$ columns.\\
{\tt simplex\_reduce} takes as argument {\tt A,b,c} et
returns  {\tt max(c.x)}, the augmented solution of {\tt x}
(augmented since the algorithm works by adding rows($A$) auxiliary
variables) and the reduced matrix.\\ 
{\bf Example}\\
Find \[ \mbox{max}(X+2Y)  \mbox{ where }
\left\{
\begin{array}{rcl}
(X,Y) & \geq & 0 \\
-3X +2Y  & \leq & 3\\
X +Y  & \leq & 4
\end{array} 
\right.
\]
Input :
\begin{center}{\tt simplex\_reduce([[-3,2],[1,1]],[3,4],[1,2])}\end{center}
Output :
\begin{center}{\tt 7,[1,3,0,0],[[0,1,1/5,3/5,3],[1,0,(-1)/5,2/5,1], [0,0,1/5,8/5,7]]}\end{center}
Which means that the maximum of {\tt X+2Y} under these conditions
is {\tt 7}, it is obtained for {\tt X=1,Y=3} 
because {\tt [1,3,0,0]} is the augmented solution and the reduced matrix is :\\
{\tt [[0,1,1/5,3/5,3],[1,0,(-1)/5,2/5,1], [0,0,1/5,8/5,7]]}.

{\bf A more complicate case that reduces to the simple case}\\
With the former call of {\tt simplex\_reduce}, we have to :
\begin{itemize}
\item rewrite constraints to the form $x_k \geq 0$,
\item remove variables without constraints,
\item add variables such that all the constraints have positive components.
\end{itemize}
For example, find :
\begin{equation}\label{eq:lpexample}
\mbox{min}(2x+y-z+4)  \quad \mbox{ where }
\left\{
\begin{array}{rcl}
x & \leq & 1 \\
y & \geq & 2 \\
x+3y-z & = & 2 \\
2x-y+z & \leq & 8\\
-x+y & \leq & 5
\end{array} 
\right.
\end{equation}
Let $x=1-X$, $y=Y+2$, $z=5-X+3Y$
the problem is equivalent to finding the minimum of
$(-2X+Y-(5-X+3Y)+8)$ 
where :
\[ 
\left\{
\begin{array}{rcl}
X & \geq & 0 \\
Y & \geq & 0 \\
2(1-X)-(Y+2)+ 5-X+3Y & \leq & 8\\
-(1-X) +(Y+2)  & \leq & 5
\end{array} 
\right.
\]
or to find the minimum of~:
\[ (-X-2Y+3) \quad \mbox{ where } 
\left\{
\begin{array}{rcl}
X & \geq & 0 \\
Y & \geq & 0 \\
-3X+2Y & \leq & 3\\
X +Y  & \leq & 4
\end{array} 
\right.
\]
i.e. to find the maximum of $-(-X-2Y+3)=X+2Y-3$
under the same conditions, hence it is the same problem as 
to find the maximum of $X+2Y$ seen before. We found {\tt 7}, 
hence, the result here is {\tt 7-3=4}.

{\bf The general case}\\
A linear programming problem may not in general be directly
reduced like above to the simple case. The reason is that
a starting vertex must be found before applying the simplex
algorithm. Therefore,
{\tt simplex\_reduce} may be called by specifying this starting
vertex, in that case, all the arguments including the starting
vertex are grouped in a single matrix. 

We first illustrate this kind
of call in the simple case where the starting point does not
require solving an auxiliary problem.
If {\tt A} has $p$ rows and $n$ columns and if we define :
\begin{center}
{\tt B:=augment(A,idn(p));} {\tt C:=border(B,b);} \\
{\tt d:=append(-c,0\$(p+1));} {\tt D:=augment(C,[d]);}
\end{center}
{\tt simplex\_reduce} may be called with {\tt D} as single argument.\\
For the previous example, input :
\begin{center}{\tt A:=[[-3,2],[1,1]];B:=augment(A,idn(2)); C:=border(B,[3,4]);
D:=augment(C,[[-1,-2,0,0,0]])}\end{center}
Here 
{\tt C=[[-3,2,1,0,3],[1,1,0,1,4]]}\\
and {\tt D=[[-3,2,1,0,3],[1,1,0,1,4],[-1,-2,0,0,0]]}\\
Input :
\begin{center}{\tt simplex\_reduce(D)}\end{center}
Output is the same result as before.

{\bf Back to the general case.}\\
The standard form of a linear programming problem is similar
to the simplest case above, but with $Ax=b$ (instead of $Ax\leq b$)
under the conditions $x\geq 0$. We may further assume that $b\geq 0$
(if not, one can change the sign of the corresponding line).
\begin{itemize}
\item The first problem is to find an $x$ in the $Ax=b, x\geq 0$ domain.
Let $m$ be the number of lines of $A$. Add artificial variables
$y_1,...,y_m$ and maximize
$-\sum y_i$ under the conditions $Ax=b, x \geq 0, y \geq 0$ 
starting with initial value $0$ for $x$ variables
and $y=b$
(to solve this with {\tt Xcas}, call \verb|simplex_reduce| with
a single matrix argument obtained by augmenting $A$ by the
identity, $b$ unchanged and an artificial
$c$ with 0 under $A$ and 1 under the identity).
If the maximum exists and is 0, the identity submatrix above the last
column corresponds to an $x$ solution, we may forget the artificial
variables (they are 0 if the maximum is 0).
\item Now we make a second call to \verb|simplex_reduce|
with the original $c$ and the value of $x$ we found in the domain.
\item
Example~: find the minimum of $2x+3y-z+t$ with 
$x,y,z,t\geq 0$ and~:
\[ \left\{ \begin{array}{rcl}
-x-y+t&=&1\\
y-z+t&=&3
\end{array}
\right. \]
This is equivalent to find the opposite of the maximum of $-(2x+3y-z+t)$.
Let us add two artificial variables $y_1$ and $y_2$,
\begin{verbatim}
simplex_reduce([[-1,-1,0,1,1,0,1],
[0,1,-1,1,0,1,3],
[0,0,0,0,1,1,0]])
\end{verbatim}
Output: optimum=0, artificial variables=0, and the matrix
\[
\left(\begin{array}{ccccccc}
-1/2 & 0 & -1/2 & 1 & 1/2 & 1/2 & 2 \\
1/2 & 1 & -1/2 & 0 & -1/2 & 1/2 & 1 \\
0 & 0 & 0 & 0 & 1 & 1 & 0
\end{array}\right) 
\]
Columns 2 and 4 are the columns of the identity (in lines 1 and 2).
Hence $x=(0,1,0,2)$ is an initial point in the domain.
We are reduced to solve the initial problem, after replacing the
lines of $Ax=b$ by the two first lines of the answer above,
removing the last columns corresponding to the artificial variables.
We add $c.x$ as last line
\begin{verbatim}
simplex_reduce([[-1/2,0,-1/2,1,2],
[1/2,1,-1/2,0,1],[2,3,-1,1,0]])
\end{verbatim}
Output: maximum=-5, hence the minimum of the opposite is 5,
obtained for $(0,1,0,2)$, after replacement 
$x=0$, $y=1$, $z=0$ and $t=2$.
\end{itemize}

For more details, search google for \verb|simplex algorithm|.

\subsection{Solving general linear programming problems: {\tt lpsolve}}\index{lpsolve|textbf}
Linear programming problems where a multivariate linear function should be maximized or minimized subject to linear equality or inequality constraints, as well as mixed integer problems, can be solved with the function {\tt lpsolve}. It takes (at most) four arguments, in the following order :
\begin{itemize}
\item {\tt obj} : symbolic expression representing the objective function,
\item {\tt constr} (optional) : list of linear constraints which may be equalities or inequalities or bounded expressions entered as {\tt expr=a..b},
\item {\tt bd} (optional) : sequence of expressions of type {\tt var=a..b} specifying that the variable {\tt var} is bounded with {\tt a} below and with {\tt b} above,
\item {\tt opts} (optional) : sequence of {\tt option=value} parameters for the solver, where {\tt option} may be one of {\tt assume}, {\tt lp\_maximize}, {\tt lp\_variables}, {\tt lp\_integervariables}, {\tt lp\_binaryvariables} or {\tt lp\_depthlimit}.
\end{itemize}
The return value is in the form {\tt [optimum,soln]} where {\tt optimum} is the minimum/maximum value of the objective function and {\tt soln} is a list of coordinates corresponding to the point at which the optimal value is attained. If there is no feasible solution, an empty list is returned. When the objective function is unbounded, {\tt optimum} is returned as {\tt +infinity} (for maximization problems) or {\tt -infinity} (for minimization problems), while {\tt soln} is generally meaningless.

If {\tt obj} is given as constant (for example, zero) then only a feasible point is returned as a list, if one exists. If the problem is infeasible, an empty list is returned. This may be used as a test to check whether a set of linear constraints is feasible or not.

The given objective function is minimized by default. To maximize it, include the option {\tt lp\_maximize=true} or simply {\tt lp\_maximize}.

By default, all variables are considered continuous and unrestricted in sign.

The solver combines the two-phase simplex method and the dual simplex method to find the optimal solution.

\subsubsection{Problems with continuous variables}
For example, to solve the problem specified in~\eqref{eq:lpexample}, input :
\begin{center}
{\tt constr:=[x<=1,y>=2,x+3y-z=2,3x-y+z<=8,-x+y<=5];}\\
{\tt lpsolve(2x+y-z+4,constr)}
\end{center}
Output :
\begin{center}
\tt [-4,[x=0,y=5,z=13]]
\end{center}
Therefore, the minimum value of $ f(x,y,z)=2\,x+y-z+4 $ is equal to $ -4 $ under the given constraints. The optimal value is attained at point $ (x,y,z)=(0,5,13) $.

Constraints may also take the form {\tt expr=a..b} for bounded linear expressions.

\noindent Input :
\begin{center}
{\tt lpsolve(x+2y+3z,[x+y=1..5,y+z+1=2..4,x>=0,y>=0])}
\end{center}
Output :
\begin{center}
{\tt [-2,[x=0,y=5,z=-4]]}
\end{center}

Use the {\tt assume=lp\_nonnegative} option to specify that all variables are nonnegative. That is easier than entering the nonnegativity constraints explicitly.

\noindent Input:
\begin{center}
{\tt lpsolve(-x-y,[y<=3x+1/2,y<=-5x+2],}\\
{\tt assume=lp\_nonnegative)}
\end{center}
Output:
\begin{center}
{\tt [-5/4,[x=3/16,y=17/16]]}
\end{center}

Bounds can be added separately for some variables. They are entered after the list of constraints.

\noindent Input :
\begin{center}
{\tt constr:=[5x-10y<=20,2z-3y=6,-x+3y<=3];}\\
{\tt lpsolve(-6x+4y+z,constr,x=1..20,y=0..inf)}
\end{center}
Output :
\begin{center}
{\tt [-133/2,[x=18,y=7,z=27/2]]}
\end{center}

\subsubsection{Integer programming}
Use the {\tt assume=integer} or {\tt assume=lp\_integer} option to solve integer programming problems. The function {\tt lpsolve} uses branch and bound method in such cases. The total number of investigated nodes is printed out before the function returns. To limit branch and bound tree depth, use the option : \begin{center}
{\tt lp\_depthlimit=<positive integer>}
\end{center} In that case optimality of the solution can't be guaranteed.

\noindent Input :
\begin{center}
{\tt lpsolve(-5x-7y,[7x+y<=35,-x+3y<=6],assume=integer)}
\end{center}
Output :
\begin{center}
{\tt [-41,[x=4,y=3]]}
\end{center}

Use the option {\tt assume=lp\_binary} to specify that all variables are binary, i.e.~the only allowed values are 0 and 1. Binary variables usually represent {\tt true}  and {\tt false} values, giving them a certain meaning in logical context.

\noindent Input :
\begin{center}
{\tt lpsolve(8x1+11x2+6x3+4x4,[5x1+7x2+4x3+3x4<=14],}\\
{\tt assume=lp\_binary,lp\_maximize)}
\end{center}
Output :
\begin{center}
{\tt [21,[x1=0,x2=1,x3=1,x4=1]]}
\end{center}

Options \begin{center}
{\tt lp\_integervariables=<list of integer variables>}
\end{center} and \begin{center}
{\tt lp\_binaryvariables=<list of binary variables>}
\end{center} are used for specifying mixed integer/binary programming problems. Also, the \begin{center}
{\tt lp\_variables=<list of continuous variables>}
\end{center} option may be used, which overrides integer and binary settings. For example, a linear programming problem with mostly integer variables may be specified using the option {\tt assume=integer} and specifying continuous variables with {\tt lp\_variables}, which overrides the global integer setting.

\noindent Input :
\begin{center}
{\tt lpsolve(x+3y+3z,[x+3y+2z<=7,2x+2y+z<=11],}\\
{\tt assume=lp\_nonnegative,lp\_maximize,}
{\tt lp\_integervariables=[x,z])}
\end{center}
Output:
\begin{center}
{\tt [10,[x=0,y=1/3,z=3]]}
\end{center}

Use the {\tt assume=lp\_nonnegint} or {\tt assume=nonnegint} option to get nonnegative integer values.

\noindent Input :
\begin{center}
{\tt lpsolve(2x+5y,[3x-y=1,x-y<=5],assume=nonnegint)}
\end{center}
Output :
\begin{center}
{\tt [12,[x=1,y=2]]}
\end{center}

\subsubsection{Entering linear programs in matrix form}
The function {\tt lpsolve} supports entering linear programming problems in matrix form, which is convenient for problems with relatively large number of variables and/or constraints.

To enter a problem in matrix form, the parameter {\tt obj} must be a vector of coefficients $ \mathbf{c} $ and {\tt constr}, which is mandatory this case, must be a list $ [\mathbf{A},\mathbf{b},\mathbf{A}_\mathrm{eq},\mathbf{b}_\mathrm{eq}] $, where $ \mathbf{A},\mathbf{A}_\mathrm{eq} $ are matrices and $ \mathbf{b},\mathbf{b}_\mathrm{eq} $ are vectors of real numbers. Without any other parameters, this minimizes $ \mathbf{c}^T\,\mathbf{x} $ under conditions $ \mathbf{A}\,\mathbf{x}\leq\mathbf{b} $ and $ \mathbf{A}_\mathrm{eq}\,\mathbf{x}=\mathbf{b}_\mathrm{eq} $. If a problem does not contain equality constraints, parameters $ \mathbf{A}_\mathrm{eq} $ and $ \mathbf{b}_\mathrm{eq} $ may be omitted. For a problem that does not contain inequality constraints, empty lists must be passed as $ \mathbf{A} $ and $ \mathbf{b} $. Also, {\tt constr} may be an empty list.

The parameter {\tt bd} is entered as a list of two vectors $ \mathbf{b}_l $ and $ \mathbf{b}_u $ of the same length as the vector $ \mathbf{c} $ such that $ \mathbf{b}_l\leq\mathbf{x}\leq\mathbf{b}_u $. For unbounded variables use {\tt +infinity} or {\tt -infinity}.

When specifying mixed problems in matrix form, variables are entered as the corresponding indexes, which are 1-based, i.e.~the first variable has index 1, the second variable has index 2 and so on. Other options for {\tt lpsolve} are passed to a problem in matrix form in the same way as if it was in symbolic form.

\noindent Input :
\begin{center}
{\tt c:=[-2,1];A:=[[-1,1],[1,1],[-1,0],[0,-1]];}\\
{\tt b:=[3,5,0,0];lpsolve(c,[A,b])}
\end{center}
Output :
\begin{center}
{\tt [-10,[5,0]]}
\end{center}
Input :
\begin{center}
{\tt c:=[-2,5,-3];bl:=[2,3,1];bu:=[6,10,3.5];}\\
{\tt lpsolve(c,[],[bl,bu])}
\end{center}
Output :
\begin{center}
{\tt [-7.5,[6.0,3.0,3.5]]}
\end{center}
Input :
\begin{center}
{\tt c:=[4,5];Aeq:=[[-1,1.5],[-3,2]];beq:=[2,3];}\\
{\tt lpsolve(c,[[],[],Aeq,beq])}
\end{center}
Output :
\begin{center}
{\tt [5.2,[-0.2,1.2]]}
\end{center}
Input :
\begin{center}
{\tt c:=[2,-3,-5];A:=[[-5,4,-5],[2,5,7],[2,-3,4]];}\\
{\tt b:=[3,1,-2];lpsolve(c,[A,b],lp\_integervariables=[1,3])}
\end{center}
Output :
\begin{center}
{\tt [19,[1,3/4,-1]]}
\end{center}

parisse
Messages : 5734
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

Re: linear programming

Message par parisse » mer. août 23, 2017 10:15 am

lukamar a écrit :Yes, it works indeed.
I have one small syntactical improvement for lpsolve. After the line 1022 it should be inserted

Code : Tout sélectionner

vector<int> vc=findvars(cvars,vars);
for (vector<int>::const_iterator it=vc.begin();it!=vc.end();++it) {
  vartype[*it]=0;
}
I hope I have added it at the right place, it's probably safer that you resend the file (zipped).
Also, I noticed that lpsolve returns +infinity for unbounded cases of minimization problems. The correct value should be -infinity. That has troubled me before but I thought that assigning minus_inf to soln[0] (see line 1048) is OK. How to fix it? (I could not test it on my computer because of the mysterious segfaulting when using plus_inf, minus_inf, zero etc., which I already reported and hadn't been able to overcome since.)
Can you give an example of lpsolve command that returns a wrong +inf?

lukamar
Messages : 331
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

Re: linear programming

Message par lukamar » mer. août 23, 2017 10:41 am

parisse a écrit :Can you give an example of lpsolve command that returns a wrong +inf?
lpsolve(x+y,[x+y<=0])
Pièces jointes
lpsolve-fixed2.zip
(13.3 Kio) Téléchargé 231 fois

parisse
Messages : 5734
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

Re: linear programming

Message par parisse » mer. août 23, 2017 2:46 pm

Replacing the minus_inf line by soln._VECTptr->front()=minus_inf; fixes the problem. There is probably something wrong with my operator []...

lukamar
Messages : 331
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

Re: linear programming

Message par lukamar » lun. sept. 25, 2017 8:54 pm

Hello,

I have improved the lpsolve command. Many parts of code were rewritten or simplified. Also, there are some new features:
1. When solving (mixed) integer LP using branch&bound, suitable Gomory mixed integer (GMI) cuts are added to each subproblem, which usually makes computation significantly faster because less subproblems have to be investigated.
2. There is new option, 'lp_integertolerance=<positive number>'. If specified, a number x such that |x-round(x)|<integertolerance is considered to be an integer. This characterization is used when solving ILPs.
3. Another new option is 'lp_method=lp_simplex or lp_interiorpoint'. When lp_method is set to lp_interiorpoint (it is set to lp_simplex by default), the solver uses primal-dual affine scaling method which operates in floating-point arithmetic. The main function is 'primal_dual_affine_scaling'. This implementation is rather simple but still requires a relatively small number of iterations to converge. But, so far (I tested it on small and small-to-moderate size problems) it's no faster than simplex method and probably needs extra simplifications and heuristics. The computational core is solving the system of linear equations for dp (lines 456-480 in lpsolve.cc). Because its matrix is symmetric, I optimized the process of solving by using Cholesky decomposition (I used the piece of code you wrote in vecteur.cc to compute it). The option 'lp_initialpoint=<list of coordinates>' may be used to specify initial point for the algorithm. If not specified, an initial point is automatically generated.
Unfortunately my interior point implementation does not work well in the combination with branch&bound, there is more work to be done in that respect.

Please, Bernard, move enum entries from lpsolve.h to the proper place (in dispatch.h, I beileve). The reserved words beginning with lp_ that should correspond to these entries are commented on the same lines.
Pièces jointes
lpsolve-fix4.zip
(17.55 Kio) Téléchargé 235 fois

Répondre