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