share/simplex/minimize_sx.mac

   1 /*****************************************************************************
   2  *                                                                           *
   3  * ************************************************************************* *
   4  * ***                                                                   *** *
   5  * ***                         ~*~ SIMPLEX ~*~                           *** *
   6  * ***                                                                   *** *
   7  * ***               A simple implementation of the simplex              *** *
   8  * ***             algorithm for Linear Programming for Maxima.          *** *
   9  * ***                                                                   *** *
  10  * ***   This file provides functions minimize_lp and maximize_lp. This  *** *
  11  * ***   file is part of the simplex package for Maxima.                 *** *
  12  * ***                                                                   *** *
  13  * ***                                                                   *** *
  14  * ***   Copyright:  Andrej Vodopivec <andrejv@users.sourceforge.net>    *** *
  15  * ***   Version:    1.01                                                *** *
  16  * ***   License:    GPL                                                 *** *
  17  * ***                                                                   *** *
  18  * ************************************************************************* *
  19  *                                                                           *
  20  * Demo                                                                      *
  21  * =====                                                                     *
  22  *                                                                           *
  23  * 1) We want to minimize x with constraints y>=x-1, y>=-x-1, y<=x+1, y<=1-x *
  24  *    and y=x/2                                                              *
  25  *                                                                           *
  26  * Solution:                                                                 *
  27  *                                                                           *
  28  * load("simplex");                                                          *
  29  * minimize_lp(x, [y>=x-1, y>=-x-1, y<=x+1, y<=1-x, y=x/2]);                 *
  30  * => [-2/3, [x=-2/3, y=-1/3]]                                               *
  31  *                                                                           *
  32  *                                                                           *
  33  * 2) If any variable is known to be positive, you should add an optional    *
  34  *    argument to minimize_lp/maximize_lp.                                   *
  35  *    We want to maximize x+y subject to x>=0, y>=0, y<=-x/2+3, y<=-x+4.     *
  36  *                                                                           *
  37  * Solution:                                                                 *
  38  *                                                                           *
  39  * maximize_lp(x+y, [y<=-x/2+3, y<=-x+4], [x, y])                            *
  40  * => [4, [x = 2, y = 2]]                                                    *
  41  *                                                                           *
  42  *****************************************************************************/
  43
  44 define_variable(nonnegative_lp,   false, boolean,
  45                "Assume all variables are non-negative!")$
  46 alias(nonegative_lp, nonnegative_lp)$
  47 define_variable(return_input_lp, false, boolean,
  48                "Return the input to linear program, not solution!")$
  49
  50 /*****************************************************************************
  51  *                                                                           *
  52  * The minimize_lp function.                                                 *
  53  *                                                                           *
  54  *****************************************************************************/
  55
  56 minimize_lp(target, constraints, [pozitive]) := block(
  57   [var : [], A, b, c, inequalities:0, count, eq, t:expand(target), sol, s,
  58    nonpozitive : 0, j, tmpvar, keepfloat:true],
  59
  60 /*******************************************************************
  61  * Get the list of variables                                       *
  62  *******************************************************************/
  63
  64   for e in constraints do (
  65     tmpvar : listofvars(e),
  66     for v in tmpvar do
  67       if not(member(v, var)) then var : cons(v, var),
  68     if op(e)#"=" then inequalities : inequalities+1
  69   ),
  70
  71   if length(pozitive)>0 then (
  72     if listp(pozitive[1]) then pozitive : pozitive[1]
  73     else if pozitive[1]='all then pozitive : copylist(var)
  74   )
  75   else if nonnegative_lp=true then
  76     pozitive : copylist(var),
  77
  78   for v in var do
  79     if not(member(v,pozitive)) then nonpozitive : nonpozitive+1,
  80
  81 /*******************************************************************
  82  * Setup A and b for linear program                                *
  83  *******************************************************************/
  84
  85   b : makelist(0, i, 1, length(constraints)),
  86   A : zeromatrix(length(constraints), length(var)+nonpozitive+inequalities),
  87   count : 0,
  88   for i:1 thru length(constraints) do (
  89     eq : lhs(part(constraints,i))-rhs(part(constraints,i)),
  90     j : 1,
  91     for v in var do (
  92       A[i,j] : ratcoeff(eq, v),
  93       if not(constantp(A[i,j])) then
  94         error("Error: constraint not linear (1).",A[i,j]),
  95       eq : subst(v=0,eq),
  96       j : j+1,
  97       if not(member(v,pozitive)) then (
  98         A[i,j] : -A[i,j-1],
  99         j : j+1
 100       )
 101     ),
 102     if op(constraints[i])="<=" or op(constraints[i])="<" then (
 103       count : count+1,
 104       A[i, length(var)+nonpozitive+count] : 1
 105     )
 106     else if op(constraints[i])=">=" or op(constraints[i])=">" then (
 107       count : count+1,
 108       A[i, length(var)+nonpozitive+count] : -1
 109     )
 110     else if op(constraints[i])#"=" then
 111       error("Error: not a proper constraint:", constraints[i]),
 112     b[i] : -eq,
 113     if not(constantp(b[i])) then
 114       error("Error: constraint not linear (2).",b[i])
 115   ),
 116
 117 /*******************************************************************
 118  * Setup c for linear program                                      *
 119  *******************************************************************/
 120
 121   c : makelist(0, i, 1, length(var)+nonpozitive+count),
 122   j : 1,
 123   for v in var do (
 124     c[j] : ratcoeff(t, v),
 125     if not(constantp(c[j])) then
 126       error("Error: cost function not linear."),
 127     t : subst(v=0,t),
 128     j : j+1,
 129     if not(member(v,pozitive)) then (
 130       c[j] : -c[j-1],
 131       j : j+1
 132     )
 133   ),
 134
 135   if not(constantp(t)) then
 136     error("Error: cost function not linear in constrained variables."),
 137
 138   if return_input_lp then return([A, b, c]),
 139
 140 /*******************************************************************
 141  * Solve the linear program                                        *
 142  *******************************************************************/
 143
 144   sol : linear_program(A, b, c),
 145
 146   if not(listp(sol)) then sol
 147   else (
 148     if sol[2]=-inf then [-inf, []]
 149     else (
 150       s : [],
 151       j : 1,
 152       for v in var do (
 153         if member(v,pozitive) then (
 154           s : append(s, [v=sol[1][j]]),
 155           j : j+1
 156         )
 157         else (
 158           s : append(s, [v=sol[1][j]-sol[1][j+1]]),
 159           j : j+2
 160         )
 161       ),
 162       [sol[2]+t, s]
 163     )
 164   )
 165 )$
 166
 167 maximize_lp(target, constraints, [pozitive]) := block(
 168   [sol : apply(minimize_lp, append([-target], [constraints], pozitive))],
 169   if not(listp(sol)) then sol
 170   else [-sol[1], sol[2]]
 171 )$
 172
 173
 174 /*****************************************************************************
 175 *                                                                           *
 176 * %functions are for debugging purposes.                                    *
 177 *                                                                           *
 178 *****************************************************************************/
 179
 180 %prepare_standard_form_lp(target,constraints,pozitive) := block([listconstvars : false,
 181   vars, slack_vars, nonpozitive_vars, nonpozitive_vars0, nonpozitive_subs, slack_vars0, nonpozitive_subs0, vars0],
 182   vars : listofvars(append([target],constraints)),
 183   slack_vars : map(lambda([e],gensym()), constraints),
 184   nonpozitive_vars : sublist(vars,lambda([e], not(member(e,pozitive)))),
 185   nonpozitive_subs : block([mgensym : lambda([x], gensym(printf(false,"~a_",x)))],
 186     map(lambda([x], x=mgensym(x)-mgensym(x)), nonpozitive_vars)),
 187   constraints : block([ltoreq : lambda([l,r], -l >= -r)],
 188     subst(["<"=ltoreq, "<="=ltoreq], constraints)),
 189   constraints : map(lambda([e,s], if op(e)="=" then e else lhs(e)-s-rhs(e)),constraints,slack_vars),
 190   slack_vars0 : block([cv:listofvars(constraints)],
 191     sublist(slack_vars, lambda([x],member(x,cv)))),
 192   nonpozitive_vars0 : listofvars(map('rhs,nonpozitive_subs)),
 193   vars0 : block([vars:vars,nonpozitive_subs1:map('lhs,nonpozitive_subs)],
 194     for v in nonpozitive_subs1 do vars:delete(v,vars), vars),
 195   [target, constraints] : subst("="="-", subst(nonpozitive_subs, [target,constraints])),
 196   [target, constraints, vars0, nonpozitive_subs, nonpozitive_vars0, slack_vars0]
 197   )$
 198 %minimize_lp(target, constraints, [pozitive]) := block([A, b, c, vars0, vars, tgt, vals, lp,
 199   nonpozitive_subs, nonpozitive_vars0, slack_vars0, mcoefmatrix],
 200   local(mcoefmatrix),
 201   mcoefmatrix(l,v) := apply('matrix,outermap(coeff,l,v)),
 202   pozitive : if pozitive=[] then [] else if listp(part(pozitive,1)) then part(pozitive,1) else if part(pozitive,1)='all or nonnegative_lp then copylist(block([listconstvars:false], listofvars(constraints))) else error("Unrecognized input."),
 203   [target, constraints, vars0, nonpozitive_subs, nonpozitive_vars0, slack_vars0] : %prepare_standard_form_lp(target, constraints, pozitive),
 204   vars : append(vars0, nonpozitive_vars0, slack_vars0),
 205   lp : linear_program(mcoefmatrix(constraints,vars), -subst(map(lambda([x],x=0),vars), constraints), first(args(mcoefmatrix([target],vars)))),
 206   if atom(lp) then lp else ([vals,tgt] : lp,
 207     block([maperror:false, mapprint:false],
 208       [tgt, append(map("=",vars0,vals),subst(map("=",nonpozitive_vars0,rest(vals,length(vars0))),nonpozitive_subs))])));
 209 %maximize_lp(target, constraints, [pozitive]) := block([tgt, vals, lp],
 210   lp : %minimize_lp(-target, constraints,first(pozitive)),
 211   if atom(lp) then lp else ([tgt,vals] : lp,
 212     [-tgt, vals]))$