Artifact 22abaf3ba24f11294fa194fa640d872fce85b3957936ff04bcb3dfc6e35d15dc:
- Executable file
r37/packages/solve/linineq.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 16869) [annotate] [blame] [check-ins using] [more...]
module linineq; % Linear inequalities and linear optimization. % Author: Herbert Melenk <melenk@zib-berlin.dbp.de> % Version 1 January 1990 % Version 1.1 February 1990 % added parameter "record=t" % Version 2 May 1991 % added Branch-and-Bound for Integer Prgramming % Version 3 Dec 1994 % added formal simplifier for MAX/MIN expressions. % Changed "inf" to "infinity". % Operator linineq_solve new. % Version 4 Jan 95 % use polytope points for the simplification of MAX/MIN % expressions. % % Solution of linear inequalities & equations with numerical % coefficients. % % Fourier(1826) /Motzkin(1936): George. B. Dantzig, % Linear Programming and Extensions. put('linineq,'psopfn, function (lambda(u); rederr "USE simplex (package linalg) instead")); global '(!*trlinineq !*trlinineqint !*prlinineq); switch trlinineq,prlinineq,trlinineqint; fluid '(linineqinterval!* linineqrecord!*); fluid '(!*ineqerr); % error code symbolic procedure linineqeval(u); % Interface for algebraic mode. begin scalar prob,equa,requa,vars,oldorder,res,u1,x,y,p,e,msg; scalar direction,rec,linineqrecord!*,r,intvars,w1,w2,op; msg := !*prlinineq or !*trlinineq; !*ineqerr :=nil; u1 := getrlist reval car u; u := cdr u; if u and eqcar(x:= reval car u,'LIST) then <<vars := getrlist x; u:=cdr u>>; while u do <<x := reval car u; u := cdr u; if eqcar(x,'EQUAL) and ((cadr x='RECORD and (rec:=t)) or (cadr x='INT and (intvars:=getrlist caddr x))) then t else <<!*ineqerr:=2; typerr(x,"illegal parameter")>>; >>; x := nil; for each u in vars do <<u := reval u; if eqcar(u,'equal) then if member(caddr u,'(min max)) then <<direction:=(cadr u. caddr u) . direction; u := cadr u>> else <<!*ineqerr:=2; rederr "illegal form in 2nd parameter">>; if smember(u,u1) and not member(u,x) then x:=u.x>>; x := vars := reversip x; while u1 do <<u := reval car u1; u1:=cdr u1; if not pairp u or not (car u memq '(geq leq equal)) then <<!*ineqerr:=2; typerr(u,"inequality")>>; op:=car u; w1:=reval cadr u; w2:=reval caddr u; if op='geq then if smemq('infinity,w2) then nil else if eqcar(w2,'max) then for each q in cdr w2 do u1:=append(u1,{{'geq,w1,q}}) else prob := (simp w1 . simp w2) . prob else if op='leq then if smemq('infinity,w2) then nil else if eqcar(w2,'min) then for each q in cdr w2 do u1:=append(u1,{{'leq,w1,q}}) else prob := (simp w2 . simp w1) . prob else if op='equal then if eqcar(w2,'!*interval!*) then u1:=append(u1,{{'geq,w1,cadr w2},{'leq,w1,caddr w2}}) else equa := (simp w1 . simp w2) . equa else <<!*ineqerr:=1; typerr(u,"inequality")>> >>; % control the linearity for each p in append(equa,prob) do <<if not domainp denr car p or not domainp denr cdr p then<<!*ineqerr:=1; rederr "unable to process nonlinear system">>; vars := linineqevaltest(numr car p, linineqevaltest(numr cdr p,vars)); >>; if msg then <<prin2 "variables:"; prin2t vars;>>; oldorder := setkorder vars; prob := for each p in prob collect (reorder numr car p./denr car p). (reorder numr cdr p./denr cdr p); equa:= for each p in equa collect (reorder numr car p./denr car p). (reorder numr cdr p./denr cdr p); % eliminate variables from equations while equa do <<e := car equa; equa := cdr equa; e := addsq(car e,negsq cdr e); if domainp numr e then <<if numr e then % nonzero constant equated to 0 <<!*ineqerr:=0; rederr "equation part inconsistent">> >> else <<u := list((x := mvar numr e). prepsq(y:=multsq(negf red numr e ./ 1, invsq(lc numr e ./ 1)))); if member(x,intvars) then % Dont eliminate integer variables - represent % equation by double inequality instead. <<x:=simp x; prob:=append({x.y,y.x},prob)>> else << prob := for each p in prob collect subsq(car p,u).subsq(cdr p,u); equa := for each p in equa collect subsq(car p,u).subsq(cdr p,u); requa := append(u,requa); if msg then <<prin2 " ";prin2 x; prin2 " eliminated by equation"; terpri()>>; vars := delete(x,vars); >> >> >>; res := if intvars then linineqint(prob,vars,msg,direction,rec,intvars) else linineq1(prob,vars,msg,direction,rec); % backsubstitution in equations; if null res then return '(list) else if res=t then res:=nil; for each e in requa do <<x := prepsq subsq(y:=simp cdr e,res); res := (car e . x) . res; if rec then <<x := prepsq y; linineqrecord!* := list(x,x) . linineqrecord!*>>; >>; setkorder oldorder; r := if rec then for each p in liqsimp!-maxmin pair(res,linineqrecord!*) collect list('LIST,list('EQUAL,caar p,cdar p),cadr p,caddr p) else for each p in res collect list('EQUAL,car p,cdr p); return 'LIST . r; end; % put('linineq_solve,'psopfn,'linineqseval); symbolic procedure linineqseval(u); 'list . reversip for each q in cdr linineqeval append(u,'((equal record t))) collect {'equal,cadr cadr q, if caddr q = cadddr q then caddr q else '!*interval!* . cddr q}; symbolic procedure linineqevaltest(f,v); % Collect the variables in standard form f and control linearity. if domainp f then v else if not(ldeg f=1) then <<!*ineqerr:=1; rederr "unable to process nonlinear system">> else if member(mvar f,v) then linineqevaltest(red f,v) else linineqevaltest(red f,mvar f.v); symbolic procedure linineq0(prob,vars,dir,rec); % Interface for symbolic mode. % Prob is a list (e1,e2,..) of algebraic expressions without % relational operators, which are interpreted as % set of inequalities ei >= 0. They are linear in the % variables vars. % Silent operation: result=nil if the system is inconsistent. begin scalar oldorder,res; linineqrecord!* := nil; oldorder := setkorder vars; prob := for each u in prob collect simp u . (nil ./ 1); res := linineq1(prob,vars,nil,dir,rec); setkorder oldorder; return res; end; symbolic procedure linineqint(prob,vars,msg,dir,rec,intvars); begin scalar x,x0,y,y0,y1,z,w,problems,best,z,z0,zbest,zf,bestr; % test integer variables and adjust order; for each x in vars do if member(x,intvars) then<<w:=x.w;intvars:=delete(x,intvars)>>; if intvars then <<!*ineqerr:=2;typerr('list.intvars,"int variables")>>; intvars:=reversip w; % select primary optimization principle. if dir then<<z:=caar dir;zf:=if cdar dir='max then 1 else -1>>; problems:=list (nil.prob); % macro loop. while problems do <<z0:=caar problems; prob:=cdar problems; problems:=cdr problems; if msg or !*trlinineqint then linineqprint2("=== next integer subproblem",prob); w:=if best and not evalgreaterp({'times,zf,z0},{'times,zf,zbest}) then nil % skip problem with suboptimal bound. else linineq1(prob,vars,msg,dir,rec); if !*trlinineqint then linineqprint3("=== subresult",w); if w and dir then <<% is better than best so far? z0:=cdr assoc(z,w); if best and evalgreaterp({'times,zf,zbest},{'times,zf,z0}) then w:=nil; >>; if w then <<% test feasability; y:=list prob; for each x in intvars do <<x0:=cdr assoc(x,w); if not fixp x0 then % branch and bound <<x:=simp x; y0:=simp{'ceiling,x0}; y1:=simp {'floor,x0}; y:= for each q in y join {(x.y0).q, (y1.x).q}; if msg or !*trlinineqint then <<writepri("branch and bound with",'first); writepri(mkquote{'list,{'geq,x:=prepsq x,prepsq y0}, {'leq,x,prepsq y1}},'last); >>; >>; >>; if cdr y then problems:=append(problems,for each q in y collect z0.q) else <<zbest:=z0; best:=w; bestr:=linineqrecord!*; if !*trlinineqint then prin2t "===> is feasable";>>; >>; % if w % without target dont need additional result. if best and null dir then problems:=nil; >>; % while problems linineqrecord!*:=bestr; return best; end; symbolic procedure linineq1(prob,vars,msg,dir,rec); % Algebraic evaluation of a set of inequalities: % prob is a list of pairs of standard quotients, % (( p1 . q1)(p2 . q2) .. (pn . qn)) % which are interpreted as inequalities: % pi >= qi ; % vars is the list of (linear) variables. % dir the direction of final optimization % rec switch; if t, the record of inequatlities is produced % Result is NIL if the system has no solution; otherwise % the solution has the form of an association list % ((v1 . val1)(v2 . val2) ... (vn . valn)), % where vi are the variables and vali are values in algebraic % form. NIL if the system has no solution. % begin scalar v,vq,lh,rh,x,y,z,prob1,prob2,prob3,prob4,nprob,sw,sol; if null vars then return linineq2(prob,msg); v := car vars; vars := cdr vars; vq := mksq(v,1); if !*trlinineq then linineqprint2(list("next variable:",v,"; initial system:"),prob); prob := linineqnormalize prob; for each p in prob do <<lh := car p; rh := cdr p; % if v appears on the lhs, isolate it if not domainp numr lh and mvar numr lh = v then <<x := invsq(lc numr lh ./ 1); sw := (numr x < 0); lh := multsq(lh,x); rh := multsq(rh,x); rh := addsq(rh,negf red numr lh ./ denr lh); if not sw then prob1 := (vq . rh) . prob1 else prob2 := (rh . vq) . prob2; >>else if domainp numr rh and domainp numr lh then prob4 := (lh . rh) . prob4 else prob3 := (lh . rh) . prob3; >>; if null prob1 and null prob2 and vars then << sol := linineq1(prob,vars,msg,dir,rec); if rec then linineqrecord!* := append(linineqrecord!*,'(((minus infinity) infinity))); return if sol then (v . 0) . sol else nil>>; if !*trlinineq then <<linineqprint2("class 1:",prob1); linineqprint2("class 2:",prob2); linineqprint2("class 3:",prob3); linineqprint2("class 4:",prob4); >>; if rec then << x := for each u in prob1 collect prepsq cdr u; y := for each u in prob2 collect prepsq car u; x := if null x then '(minus infinity) else if null cdr x then car x else 'max . x; y := if null y then ' infinity else if null cdr y then car y else 'min . y; linineqrecord!* := append(linineqrecord!*, list list(x,y)) >>; if not linineq2(prob4,msg) then return nil; nprob := append(prob3, for each x in prob1 join for each y in prob2 collect car y . cdr x); if vars then << if null (sol := linineq1(nprob,vars,msg,dir,rec)) then return nil>> else if not linineq2(nprob,msg) then return nil; % lower bound: x := if null prob1 then nil else linineqevalmax for each p in prob1 collect subsq(cdr p,sol); % upper bound: y := if null prob2 then nil else linineqevalmin for each p in prob2 collect subsq(car p,sol); if (z:=assoc(v,dir)) then z:= cdr z; if msg then <<writepri(" ",'first); writepri(mkquote if x then prepsq x else '(minus infinity),nil); writepri(" <= ",nil); writepri(mkquote v,nil); writepri(" <= ",nil); writepri(mkquote if y then prepsq y else 'infinity,nil); writepri("; ",nil)>>; linineqinterval!* := x . y; if z='min and null x or z='max and null y then <<if msg then writepri( " max/min cannot be resolved",'last); return nil>>; if not(x=y) then if z='min then y:=nil else if z='max then x:=nil; if msg then << writepri( if null x and null y then " completely free: " else if null y then " minimum: " else if null x then " maximum: " else if x=y then " zero length interval: " else " middle: ",nil); >>; if null x and null y then x := 0 else % completely free if null x then x := prepsq y else if null y then x := prepsq x else if sqlessp(y,x) then <<prin2 "system inconsistent:"; prin2 prepsq x; prin2 " not <= "; prin2t prepsq y; return nil>> else x := list('quotient,list('plus,prepsq x,prepsq y),2); x := aeval x; if msg then writepri(mkquote list('EQUAL,v,x),'last); return (v . x) . sol; end; symbolic procedure linineq2(prob,msg); % All variables are elimitated. Control, if the % remaining numerical inequalities are consistent. begin scalar rh,lh; loop: if null prob then return t; lh := caar prob; rh := cdar prob; if not domainp numr rh or not domainp numr lh then <<!*ineqerr:=1; rederr{" non numeric:", rh, lh}; >>; if sqlessp(lh,rh) then <<if msg then <<writepri("system inconsistent: ",'first); writepri(mkquote prepsq lh,nil); writepri(" not >= ",nil); writepri(mkquote prepsq rh,'last);>>; return nil>>; prob := cdr prob; goto loop; end; symbolic procedure linineqnormalize prob; % Normalize system: reform all inequalities such that they have % the canonical form % polynomial >= constant % (canonical restriction: absolute term of lhs=0, % denominator of lhs = 1). % and remove those, which have same lhs, but smaller rhs % (the latter are superfluous). begin scalar r,lh,rh,d,ab,x; for each p in prob do <<lh := car p; rh := cdr p; % arithmetic normalizaion lh := addsq(lh,negsq rh); d := denr lh; lh := numr lh; ab := lh; x := if domainp lh then 1 else lc ab; while not domainp ab do <<x := gcdf(x,lc ab);ab := red ab>>; ab := negf ab; lh := multsq(addf(lh,ab)./1,1 ./ x); rh := multsq(ab ./ 1, 1 ./ x); % removal of redundant elements x := assoc(lh,r); if null x then r:=(lh.rh) . r else if sqlessp(cdr x,rh) then rplacd(x,rh); >>; if !*trlinineq then linineqprint2("normalized and reduced:",r); return r; end; symbolic procedure linineqevalmin u; % Compute the minimum among the list u with sq's. linineqevalmin1(car u,cdr u); symbolic procedure linineqevalmin1(q,u); if null u then q else (linineqevalmin1( if x and !:minusp x then q else car u, cdr u) ) where x=numr addsq(q,negsq car u); symbolic procedure linineqevalmax u; % compute the maximum among the list u with sq's linineqevalmax1(car u,cdr u); symbolic procedure linineqevalmax1(q,u); if null u then q else (linineqevalmax1( if x and !:minusp x then car u else q, cdr u) ) where x=numr addsq(q,negsq car u); symbolic procedure sqlessp(q1,q2); (x and !:minusp x) where x=numr addsq(q1,negsq q2); symbolic procedure liqsimp!-maxmin w; liqsimp2!-maxmin liqsimp1!-maxmin w; endmodule; end;