clear //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// function Gout=UmaLineaAxLeqBx(A,B,Gin) //G=UmaLineaAxLeqBx(A,B,Gin) //A, B, et Gin donnees, Gout est une matrice dont les colonnes sont des generateurs des solutions de Ax<=Bx. //Si l'ensemble est vide, Gout est une matrice nulle (\varepsilon partout). Gin est //l'ensemble des generateus a tester. //Cette procedure provient de l'algorithme de Xavier Allamigeon 'Double sided method' 2010. // [r,n]=size(Gin); //G=%zeros(r,n);H=G; new=%F; if (full(A) > full(B)) then Gout=%zeros(r,1);//disp('Pas de solution.') elseif Gin==%zeros(r,1) then Gout=%zeros(r,1); else p=0;q=0; G=%zeros(r,1); H=%zeros(r,1);GG=%zeros(r,1);GGG=%zeros(r,1); for i=1:n x=Gin(:,i); if (full(A*x)<=full(B*x)) then q=q+1; G(:,q)=x; else p=p+1;H(:,p)=x; end end; if (p>0) then k=0; for i=1:p;for j=1:q; k=k+1; GG(:,k)=(((A*H(:,i))*G(:,j))+((B*G(:,j))*H(:,i))); end;end; j=0; G2=G(:,1:q); for i=1:k; mb=-max(plustimes(GG(:,i)));b=(GG(:,i)*mb); //v=max(plustimes(G2*(G2\b))); Este no funciona, pues G2 no es cuadrada en general if (~includespan(full(G2),full(b))) then//(v<>b) then j=j+1; GGG(:,j)=b; G2=[G(:,1:q) GGG]; new=%T; end end end if (new) then Gout=[G(:,1:q) GGG]; else Gout=G; end //Gout=G; end endfunction //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// function Gout=VariasLineasAxLeqBx(A,B,Gin) [r,n]=size(A) //trivial=%zeros(n,1); Gout=Gin; for i=1:r // QuasiVazio=%F; // i=1; // while ((i<=r) & (~QuasiVazio)) Gout=UmaLineaAxLeqBx(A(i,:),B(i,:),Gout); // if (Gout==trivial) then QuasiVazio=%T // i=i+1; // end; end endfunction //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// function [A,B]=FastCono(M) [linm,colm]=size(M); Ap=[M' %zeros(colm,linm)]; [linAp,colAp]=size(Ap); Bp=[%zeros(colm,linm) M']; Ojo=%eye(2*linm,2*linm); for i=1:linAp gleq=VariasLineasAxLeqBx(Ap(i,:),Bp(i,:),Ojo); gleq=weakbasis(full(gleq)); geq=VariasLineasAxLeqBx(Bp(i,:),Ap(i,:),gleq); Ojo=weakbasis(full(geq)); end A=geq(1:linm,:)'; B=geq(linm+1:2*linm,:)'; endfunction ///////////////////////////////////////////////////////////////////////////////////////////// function N=Twin(M) N=maxplus(-plustimes(full(M))); endfunction ///////////////////////////////////////////////////////////////////////////////////////////// function X=mpm(Z,Y) // min plus multiplication of two max-plus matrices [liny,coly]=size(Y); [linz,colz]=size(Z); Y=plustimes(Y);Z=-plustimes(Z); X=maxplus(Z)*maxplus(-Y); X=plustimes(X); X=maxplus(-X); endfunction ///////////////////////////////////////////////////////////////////////////////////////////// function PsiPlus1=PsiI(alpha,M,A,B,C,D,x) // x is the precedent psi // C and D are the matrices for the cone representation of the AB-invariant module // A and B descriB the system, x(k+1)=Ax(k)+Bu(k), A*mu is the column j of A when mu is a vector // made of %0 except at line j. // One apply this algorithm for each column of the friend F choosing the right mu // The algorithm stops when PsiPlus1 is equal to precedent PsiPlus1. This function descriBd // here is one step of the algorithm. Bta=mpm(mpm(maxplus(-plustimes(C*B)'),(D*B*x*M+D*A*M)),maxplus(-plustimes(M)')); //mpm is defined above Gamma=mpm(mpm(maxplus(-plustimes(D*B)'),(C*B*x*M+C*A*M)),maxplus(-plustimes(M)'));; PsiPlus1=maxplus(-plustimes(alpha))+maxplus(-plustimes(Bta))+maxplus(-plustimes(Gamma)); PsiPlus1=maxplus(-plustimes(PsiPlus1)); endfunction ///////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////////////// function B=Normalizep(A) //Do not work yet, well, maybe it works now, try it! [nlin,ncol]=size(A); A=Normalize(A); B=%zeros(nlin,ncol); for i=1:ncol r=Twin(mpm(%ones(1,nlin),(A(:,i)))); B(:,i)=r*A(:,i); end; endfunction /////////////////////////////////////////////////////////////////////////////////////////////// function B=Normalize(A) [nlin,ncol]=size(A); B=%zeros(nlin,ncol); for i=1:ncol r=mpm(Twin(%ones(1,nlin)*A(:,i))',%1); B(:,i)=r*A(:,i); end; endfunction /////////////////////////////////////////////////////////////////////////////////////////////// //getf('Twin.sci'); //getf('UmaLineaAxLeqBx.sci'); //getf('VariasLineasAxLeqBx.sci'); //getf('mpm3.sci'); //getf('Cono.sci'); //getf('Slim2.sci'); //getf('Normalize.sci'); //getf('Normalizep.sci'); N=50; C=[%0 %0 %0 0]; D=[0 %0 %0 0]; A=[0 -10 -10 1;0 0 -10 2;-10 0 0 3;-10 -10 0 4];A=maxplus(A); B=[0 0 0 0]';B=maxplus(B); disp(full(A),'A='); disp(full(B),'B='); disp(full(C),'C='); disp(full(D),'D='); R=A(:,4); [rho,v]=howard(A); [nlin,ncol]=size(A); [Pd,Qd]=mpkernel(D); //The kernel of D, is also a module: [Pd;Qd]. [Pdic,Qdic]=mpkernel([D;C]); //The intersection of KerC and KerD, is also a module [Pdic;Qdic]. //And A(KerC \Cap KerD) is [A*Pdic;A*Qdic], that is not always a kernel of some matrix, but here we can //take it esay ;-) disp(includespan(full([Pd;Qd]),full([A*Pdic;A*Qdic])),'Is A(KerC \Cap KerD) \subset KerD ?'); //Here start the questionable algorithm - not questionable anymore, JJ Loiseau prooved that A * CalC is also //transitive. He prooved also that if S0 is a congruence, then KerD \oplus A (S^k \Cap KerC ) is also a congruence.Mars 2025. // Well the following lines need to be verified. S0 should be KerD. I=full(%eye(nlin,ncol)); [S0g,S0d]=mpkernel(I); S0=[S0g;S0d]; //At the first step (well the one after step0), we want A*(S0 \Cap KerC) [SiCg,SiCd]=mpkernel([I;C]); S1=[[S0g;S0d] [A*SiCg;A*SiCd]]; disp(includespan(full([Pd;Qd]),S1),'Is S1 \subset KerD ?'); //So we can conclude that we can find an observer! We could also read the article ;-) ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////