//***************************************************
// TEXTE : FIABILITÉ DES SYSTÈMES
// Scopos vol 11, p. 9.
//***************************************************

// Jeu de données test (matrice A p. 16 eq. 13) :
// donnees(0.001,1,0.001,10)
function [A]=donnees(lambda1,lambda2,mu1,mu2)
//lambda1=0.001;
//lambda2=1;
//mu1=0.001;
//mu2=10;

A=zeros(4,4);
A(1,1)=-lambda1-lambda2;
A(1,2)=lambda1;
A(1,3)=lambda2;

A(2,1)=mu1;
A(2,2)=-mu1-lambda2;
A(2,4)=lambda2;

A(3,1)=mu2;
A(3,3)=-lambda1-mu2;
A(3,4)=lambda1;

A(4,2)=mu2;
A(4,3)=mu1;
A(4,4)=-mu1-mu2;

endfunction

// Question 2 : Calcul direct de la probabilité stationnaire (disponibilité),
// on rédout le système PA=0.
function [P]=Plim(A)
P=kernel(A'); P=P ./ sum(P);
endfunction
// N.B. : l'INdisponibilité d'un système de 2 composants en parallèle est le
// produit des INdisponibilités (indépendance : le système global est en panne
// sssi les 2 composants le sont). 
// On a donc une expression simple pour P en utilisant le calcul fait en cours
// pour une seule machine (processus markovien de saut à 2 états).
// Pour une machine D= mu / (lambda+mu), d'où pour les 2 machines en parallèle
// 1-P = lambda1 / (lambda1+mu1) * lambda2 / (lambda2+mu2).

// Si D1=Plim(A1) et D2=Plim(A2) sont les disponibilités obtenues pour les
// deux systèmes, celle du système série est donnée par D = D1 * D2, 
// (indépendance : le système global est en marche sssi les 2 sous-sytèmes
// le sont).

// Calcul de exp(tA) pour t fixé et A une matrice carrée
// par scaling/squaring. 
function [E]=expss(A,t)
N=10;
a=0.35;
// N.B. : pour toute norme matricielle qui vérifie 
// || A * B || <= ||A|| * ||B||, (c'est le cas de la norme
// ||A|| = max_i (\sum_j |a_ij|) appelée norm(A,'inf') en scilab)
// si  t ||A|| < a, l'erreur de troncature sur exp(tA) 
// est majorée par   a^(N+1) / (N+1)! / (1 - a/N+2), 
// soit avec a=0.35 et N=10 :
// eps = 0.35^11 / prod(1:11) / (1- 0.35 /12) < 2.5 10^{-13}
//
// Scaling : si t*||A|| > a, on divise t par une puissance de 2 telle que
// t /2^n * ||A|| < a.
if t * norm(A,'inf') > a,
   n= ceil(log2(t * norm(A,'inf') / a));
else
   n=0;
end
//
// On calcule ensuite exp(B) avec 
B= t .* A ./ (2^n);
E=eye(B)+B;
puissB=B;
factorielle=1;
for k=2:N
	puissB=puissB*B;
	factorielle=factorielle*k;
	E=E+puissB/factorielle;
end;
// E contient maintenant la série exponentielle exp(B) tronquée 
// à l'ordre N.
//
// Squaring : calcul de exp(tA) = exp(2^n B) = [exp(B)]^(2^n)
// Évolution de l'erreur lors des passages au carré successifs :
// Si ||P-Q|| < eps, et si P et Q commutent (c'est le cas ici
// car P=exp(B) et Q= sa série de Taylor tronquée), alors
// P^2 -Q^2 = (P-Q)(P+Q), d'où, avec une norme vérifiant
// || A * B || <= ||A|| * ||B||,
// ||P^2 -Q^2|| est de l'ordre de 2 ||P|| * eps, ici la norme
// de P est 1, celle de ses puissances aussi, donc au bout de
// n élévations au carré l'erreur est multipliée par 2^n,
// soit 1024 pour n=10 par exemple (perte de 3 décimales pour n=10,
// de 6 décimales pour n=20).
if n>0,
  for i=1:n,
    E=E*E;
  end;
end
endfunction

// Avec la matrice 4x4 proposée ||A||=norm(A,'inf')=20
// donc jusqu'à t=10000, ||tA|| <= 200000, ce qui donne n=20
// et une erreur maxi < 10^(-6) sur exp(tA) pour t<=10000.
// La méthode permet donc de tracer le graphe de la page 19.

// Question 3.
// Tracé des graphes de la disponibilité D(t) et de la fiabilité R(t)
// du système d'embouteillage (système série) :
// Disp(t) = exp(t*A_R)(1,4)  * exp(t*A_M)(1,4) ; 
// Fiab(t) = exp(t*A0_R)(1,4) * exp(t*A0_M)(1,4).
function graphe(lambda1,lambda2,mu1,mu2)
// matrice M pour les 2 machines identiques
M=donnees(lambda1,lambda1,mu1,mu1);
// M avec état 4 absorbant :
M0=M; M0(4,:)=[0 0 0 0];
// matrice R pour les 2 réservoirs identiques
R=donnees(lambda2,lambda2,mu2,mu2);
// A avec état 4 absorbant :
R0=R; R0(4,:)=[0 0 0 0];
// T = vecteur colonne des temps
T=[0.1 0.2 0.4 0.8 1 2 4 8 10 20 40 80 100 200 400 800 1000 2000 4000 8000]';
Disp=ones(T);Fiab=ones(T);
k=1;
for i=-1:3,
  E=expss(R,10^i); E0=expss(R0,10^i);
  F=expss(M,10^i); F0=expss(M0,10^i);
  Disp(k)=(1-E(1,4))*(1-F(1,4));
  // On prend la partie positive pour les erreurs d'arrondi...
  Fiab(k)=max((1-E0(1,4))*(1-F0(1,4)), 0);
  for j=1:3,
    E=E*E; E0=E0*E0; F=F*F; F0=F0*F0;
    Disp(k+j)=(1-E(1,4))*(1-F(1,4));
    Fiab(k+j)=max((1-E0(1,4))*(1-F0(1,4)), 0);
  end
  k=k+4;
end
P=Plim(R);Q=Plim(M);
// La disponibilité du système global en régime stationnaire est
DispLim=(1-P(4))*(1-Q(4));
// Tracé des 2 courbes Disp(t) Fiab(t) et de l'asymptote DispLim
xbasc();
xselect();
// Tracé normal (échelle linéaire sur Ox et Oy)
//plot2d(T,[DispLim*ones(T),Disp,Fiab], style=[2,13,5]);
// Avec échelle logarithmique sur Ox :
plot2d(T,[DispLim*ones(T),Disp,Fiab], style=[2,13,5],logflag="ln");
endfunction

// Question subsidiaire.
// Comparer expss(A) à avec la fonction scilab expm(A) n'est pas probant, 
// car scilab utilise aussi la méthode scaling/squaring, mais ça permet au
// moins de vérifier que la fonction expss ne donne pas des résultats farfelus.
// expm(t*A) : calculer
// expm(A)-expss(A,1)
// expm(100*A)-expss(A,100)
// expm(10000*A)-expss(A,10000)

// Autre comparaison : diagonalisation sous scilab.
// spec(A) = [0, -11.002, -11, -0.002]
// donc 4 valeurs propres distinctes, A est diagonalisable.
// On a kernel(A)= [ 1; 1; 1; 1]
// kernel(A+11.002 *eye(A)) = 
//   [ 0.0703598; -0.0703598; -0.7035975; 0.7035975]
// kernel(A+11 *eye(A)) =
//   [-0.0703598; -0.0703598;  0.7035975; 0.7035975]
// kernel(A+ 0.002 *eye(A)) =
//   [ -1; 1; -1; 1]
// D'où la matrice de passage P =
// !   1.    0.0703598  - 0.0703598  - 1. !
// !   1.  - 0.0703598  - 0.0703598    1. !
// !   1.  - 0.7035975    0.7035975  - 1. !
// !   1.    0.7035975    0.7035975    1. !
//  D  =
// !   0.    0.        0.     0.    !
// !   0.  - 11.002    0.     0.    !
// !   0.    0.      - 11.    0.    !
// !   0.    0.        0.   - 0.002 !
// On vérifie que P * D * inv(P) - A est (presque) la matrice nulle
// (à 10^{-6} près...). La méthode numérique de diagonalisation
// n'est pas exempte d'approximations non plus...
// Il reste à comparer expss(A,t) à
// P * (exp(t*D) -ones(D) + eye(D)) * inv(P),
// les écarts sont inférieurs à 10^{-7}.

// Autre exemple critique donné par Moler & van Loan (Siam Review, vol. 20,
// 1978, p. 801-836) : A =[-49 24 ; -64 31];
// A est diagonalisable (v.p. -1 et -17)
// e^A = [1 3;2 4] * [e^(-1) 0; 0 e^(-17)] * inv([1 3; 2 4])
//     = [- 0.7357588    0.5518191; - 1.4715176    1.1036382]