Calcul formel : Mode d'emploi - Exemples en Maple
Claude Gomez, Bruno Salvy, Paul Zimmermann
Masson, 1995
Chapitre V, section 2.5, exercice 4, page 137.
Philippe.Dumas@inria.fr
http://algo.inria.fr/dumas/Maple/
|
|
La difficulté de l'exercice n'est pas dans le fait de récupérer les coefficients à partir des expressions. Le problème est de vérifier la correction des données et leur format, ou plutôt de se ramener à un format considéré comme standard. Cette tâche étant effectuée il suffit essentiellement d'employer coeff.
L'exercice fait référence au package simplex. Il n'est pas inintéressant de jeter un oeil sur les procédures de package. On voit en particulier comment la correction des arguments est testée.
> ?simplex
> interface(verboseproc=2):
La procédure simplex[minimize] ne fait qu'appeler simplex[maximize], ce qui n'est guère étonnant.
> print(simplex[minimize]);
Voyons simplex[maximize]. Notez le readlib(`simplex/init`)(), la vérification de type et l'appel à simplex[feasible]. Dans la suite de la page, les points remarquables du code sont marqués par un bouton .
> print(simplex[maximize]);
Nous regardons aussi simplex[feasible]. Nous voyons une fois de plus le readlib(`simplex/init`)(). Nous constatons que les inégalités strictes ne sont pas traitées. Les contraintes sont ensuite traitées : les contraintes de la forme 0 R 0 où R est l'égalité ou une inégalité large sont éliminées ; celles de la forme c R 0 signifient que le système ne possède pas de solution ; enfin on teste si celles qui restent sont de type linéaire.
> print(simplex[feasible]);
Regardons ce que comporte simplex/init. Nous voyons que le type simplex/linear est défini par une procédure simplex/type/linear.
> readlib(`simplex/init`)();
> print(`simplex/init`);
La procédure simplex/type/linear vérifie qu'une expression est de type linéaire, ou plus précisément que l'on a affaire à une expression ou à une relation entre expressions de type affine (il peut y avoir un terme constant). On notera dans le readlib les quotes (accent aigu ou apostrophe) et les backquotes (accent grave). Les backquotes sont là pour que le nom ne soit pas vu comme une fraction à cause des barres de division et les quotes évitent le problème qui surgirait si l'utilisateur avait affecté une valeur à ce nom.
> eval(`type/simplex/linear`);
Il ne nous reste plus qu'à écrire une procédure qui réponde à la question posée en nous inspirant des exemples précédents. La première phase (#phase1) consiste à vérifier les données. On notera que nous avons accepté les inégalités strictes, contrairement à ce qui est fait dans simplex[feasible]. Nous aurions pu permettre qu'une expression expr s'interprète comme l'égalité expr=0, mais cette possibilité n'est pas prise en compte dans l'ensemble du package ; nous ne l'avons donc pas permise. Enfin l'objectif est censé être linéaire, d'après l'énoncé, et nous avons seulement vérifié qu'il est affine. La seconde phase (#phase2) est un simple remplissage de matrices. On notera que nous renvoyons une liste constituée des variables qui figurent dans le critère (l'objectif) et dans les contraintes. Cette liste permet de savoir dans quel ordre sont rangées les variables. Sans cette liste, les matrices n'ont pas de sens. Cette liste a été obtenue grâce à la procédure indets.
> formalize:=proc(obj,constraints)
global `type/simplex/linear`;
local cnsts,c,x,X,Z,z,n,m,A,b,i,j;
`type/simplex/linear`:=readlib('`simplex/type/linear`');
if not type(obj,{constant, `simplex/linear`}) then #phase1
ERROR("non-linear objective function",obj,constraints) fi;
if not type(constraints,{list, set}(relation)) then
ERROR("non valid constraints",obj,constraints) fi;
cnsts:=convert(constraints,list);
cnsts:=map(proc(eq) op(1,eq)-op(2,eq) end,cnsts);
cnsts := select(proc (x) not evalb(x=0) end,cnsts);
if select(proc (x) type(x,constant) end,cnsts) <> [] then
ERROR("incompatible constraints",obj,constraints) fi;
for c in cnsts do
if not type(c,'`simplex/linear`') then
ERROR("non-linear constraint",c) fi od;
X:=indets(cnsts);
Z:=indets(obj);
for z in Z do
if not member(z,X) then
ERROR("some variables in the objective function are
not constrained",obj,constraints)
fi
od;
X:=convert(X,list);
n:=nops(cnsts); #phase2
m:=nops(X);
A:=array(1..n,1..m);
b:=array(1..n);
c:=array(1..m);
for i to n do
for j to m do
A[i,j]:=coeff(cnsts[i],X[j])
od;
b[i]:=add(A[i,j]*X[j],j=1..m)-cnsts[i]
od;
for j to m do
c[j]:=coeff(obj,X[j])
od;
X,eval(A),eval(b),eval(c)
end:
Bien sûr, en V.3 il faut remplacer le add par un convert/+. D'autre part, en V.3 et V.4 les messages d'erreur sont des noms, alors qu'en V.5 on peut employer des noms ou des chaînes de caractères, puisque les deux concepts existent. Il semble plus normal d'utiliser des chaînes de caractères, c'est pourquoi les messages sont délimités par des guillemets (") et non par des backquotes (`).
> constr:={-3*x[1]+2*x[2]=-exp(1)*x[3]+2,
-x[1]+sqrt(2)*x[2]+x[4]=4,
x[1]+x[2]+x[5]=5,
2*(x[1]+x[2])+Pi*x[3]=0};
> goal:=-x[1]-2*x[2];
> formalize(goal,constr);