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
44define_variable(nonnegative_lp,   false, boolean,
45               "Assume all variables are non-negative!")$
46alias(nonegative_lp, nonnegative_lp)$
47define_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
56minimize_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
167maximize_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]))$
213