# Objectif

Dans ce court tutoriel, nous verrons brièvement les opérations de base sur les matrices et les vecteurs, ainsi que l'accès à leurs composantes.
Le propos sera étendu en attaquant des problèmes d'algèbre linéaire simple.


```julia
using LinearAlgebra
```

# 1. Opérations générales

### Caractère unicode et LaTeX

Les caractères unicodes sont admissibles en Julia. Pour les utiliser, il suffit de faire faire la commande LaTex du caractère voulu suivi de la touche Tab. Par exemple, \alpha imprime le caractère $\alpha$.

### Scalaires et constants

Un certain nombre de constantes utiles
* $pi$ ou $\pi$ $\rightarrow$ la valeur de "$pi$"
* $e$ $\rightarrow$ la valeur de "$e$"
* $inf$ $\rightarrow$ infini
* $NaN$ $\rightarrow$ indéterminé (Not a Number)
* $eps()$ $\rightarrow$ précision machine

### Opérateurs arithmétiques de base

* $+$ $\rightarrow$ addition
* $-$ $\rightarrow$ soustraction
* $\ast$ $\rightarrow$ multiplication
* $/$ $\rightarrow$ division
* $\backslash$ $\rightarrow$ division à gauche
* $^{\wedge}$ $\rightarrow$ puissance

### Opérateurs relationnels

* $==$ $\rightarrow$ test d'égalité
* $!=$ ou $\ne$ $\rightarrow$ test différence
* $<$ $\rightarrow$ test d'infériorité
* $>$ $\rightarrow$ test de supériorité
* $<=$ ou $\le$ $\rightarrow$ test d'infériorité ou égalité
* $>=$ ou $\ge$ $\rightarrow$ test de supériorité
* false $\rightarrow$ faux (logique)
* true $\rightarrow$ vrai (logique)

### Opérateurs logiques

* ! $\rightarrow$ négation logique
* expression1 & expression2 $\rightarrow$ et logique
* expression1 && expression2 $\rightarrow$ et logique
* expression1 $\parallel$ expression2  $\rightarrow$ ou logique

# 2. Matrices

### Déclarations

Crée une matrice de 2 lignes et 3 colonnes contenant les éléments 1 2 3 sur la première ligne et 4 5 6 sur la deuxième;


```julia
a = [1 2 3; 4 5 6]
```




    2×3 Matrix{Int64}:
     1  2  3
     4  5  6



Crée un vecteur colonne valant 1 2 3


```julia
b=[1 ;2 ;3]
a=[1,2,3] # opération identique
```




    3-element Vector{Int64}:
     1
     2
     3



L'opération de transposition conjuguée est réalisée avec l'opérateur ':


```julia
b=a'
```




    1×3 adjoint(::Vector{Int64}) with eltype Int64:
     1  2  3



Pour obtenir la transposée sans la congugaison, on emploiera la fonction `transpose`.


```julia
b = transpose(a)
```




    1×3 transpose(::Vector{Int64}) with eltype Int64:
     1  2  3



Il est possible de concaténer des objets avec la fonction *vcat*.

L'exemple suivant crée une matrice dont la première ligne vaut $a$ et la deuxième ligne vaut [ 1 2 3 ].


```julia
vcat(a',[1 2 3])
[a'; [1 2 3]]# opération identique
```




    2×3 Matrix{Int64}:
     1  2  3
     1  2  3



Crée un vecteur valant 1 2 3 4 5


```julia
a=1:5
a=range(1,stop=5) # opération identique
```




    1:5




```julia
println(typeof(a))
println(a) 
```

    UnitRange{Int64}
    1:5



```julia
b=collect(a)
println(typeof(b))
println(b)
```

    Vector{Int64}
    [1, 2, 3, 4, 5]


Il est possible d'utiliser un autre incrément que 1. La syntaxe est début:incrément:fin. Par exemple, la ligne suivante crée un vecteur valant [ 1 3 5 7 9 ]'


```julia
a=collect(1:2:10)
```




    5-element Vector{Int64}:
     1
     3
     5
     7
     9



Crée un vecteur valant 1 3 4 5


```julia
a=[1 collect(3:5)']
```




    1×4 adjoint(::Vector{Int64}) with eltype Int64:
     1  3  4  5



Crée une matrice de 0 de taille $3 \times 3$;


```julia
a=zeros(3, 3)
```




    3×3 Matrix{Float64}:
     0.0  0.0  0.0
     0.0  0.0  0.0
     0.0  0.0  0.0



La commande zeros est très souvent utilisée pour déclarer une matrice, vu qu'en Julia, les variables ne sont pas déclarées.

Crée une matrice de 1 de taille $2 \times 4$.


```julia
a=ones(2,4)
```




    2×4 Matrix{Float64}:
     1.0  1.0  1.0  1.0
     1.0  1.0  1.0  1.0



### Accès aux éléments et sous-matrices

Soit la matrice


```julia
A=[1 2 3;
   4 5 6]
```




    2×3 Matrix{Int64}:
     1  2  3
     4  5  6



et le vecteur


```julia
b=[1 2 3]
```




    1×3 Matrix{Int64}:
     1  2  3



L'accès à l'élément $(i,j)$ de $A$ s'obtient simplement avec $A[i,j]$.
On peut accéder à une sous-matrice en entrant l'indice de début et l'indice de fin, séparés de :.
Omettre les indices revient à sélectionner tous les indices, par exemple


```julia
A[:,2:3]
```




    2×2 Matrix{Int64}:
     2  3
     5  6



Par exemple, sélectionnons les colonnes 1 et 3 de $A$


```julia
B = [ A[:,1] A[:,3] ]
```




    2×2 Matrix{Int64}:
     1  3
     4  6



$A[:]$ crée un vecteur colonne formé par tous les éléments de $A$.

En général, si $mat$ est une matrice, alors
* $mat[i,j]$ $\rightarrow$ l'élément (i,j) de mat
* $mat[i,:]$ $\rightarrow$ la ligne i de mat
* $mat[i:j,:]$ $\rightarrow$ les lignes i à j de mat 
* $mat[:,j]$ $\rightarrow$ la colonne j de mat
* $mat[:,k:m]$ $\rightarrow$ les colonne k à m de mat
* $mat[i:j,k:m]$ $\rightarrow$ les éléments se trouvent dans laes lignes i à j et dans les colonnes k à m de mat

### Opérations de base

* A': transposition normale si $A$ est une matrice réelles et transposition conjuguée si $A$ est une matrice complexe  (conjugué de la transposée);
* C=B^-1 ou C=inv(B): inverse de B;
* C*B: produit matriciel;
* C.*B: produit terme à terme;
* C+B: somme;
* C^2: $C.*C$ (puissance matricielle);
* C.^B: Puissance terme à terme;
* C./B: division terme à terme;
* B\C: équivalent à B^-1*C;
* B/C: Equivalent à (C'\B')'.

Note: on peut remplacer une des opérandes par un scalaire. Par exemple B+2 ajoute 2 à chacun des termes de $B$.

### Fonctions usuelles sur les matrices

* size(B): vecteur correspond aux dimensions de $B$;
* size(B,1): vecteur correspond à la taille de la première dimension de $B$;
* length(v): taille de $v$;
* length(B): somme entre le nombre de lignes et le nombre de colonnes (length($B$) = sum(size($B))$);
* x=diag(B): vecteur colonne contenant la diagonale de $B$;
* diag(x): matrice diagonale dont la diagonale correspond à $x$;
* det(B): déterminant;
* norm(B): norme;
* rank(B): rang;
* trace(B): trace;
* minimum(x): minimum des composantes de $x$;
* maximum(x): maximum des composantes de $x$;

#### Exemple d'opérations de bases
On veut résoudre (pour x) l'équation Ax=b


```julia
A = [1 2;
   3 4]
b = [5;6]
x=inv(A)*b #donne la solution désiré, mais FORTEMENT déconseillé
x=A\b # analytiquement équivalent, numériquement supérieur
```




    2-element Vector{Float64}:
     -3.9999999999999987
      4.499999999999999



On veut multiplier, terme à terme, les éléments de $A$ par eux même. 


```julia
B = A .* A
B = A .^ 2 # opération identique
```




    2×2 Matrix{Int64}:
     1   4
     9  16



Il faut bien comprendre que c'est très différent de calculer $A^2 = A*A$ où * est la multiplication matricielle


```julia
C = A * A # C ≠ B
C = A ^ 2 # opération identique et toujours C ≠ B
```




    2×2 Matrix{Int64}:
      7  10
     15  22



La plupart de ces fonctions nécessitent l'utilisation de la librairie LinearAlgebra, installée par défaut avec Julia. Pour ce faire, nous entrerons


```julia
using LinearAlgebra
```

#### ll est possible d'obtenir de l'aide sur une fonction avec la commande ?, suivi du nom de la fonction.

Par exemple


```julia
?rank
```

    search: [0m[1mr[22m[0m[1ma[22m[0m[1mn[22m[0m[1mk[22m [0m[1mR[22m[0m[1ma[22m[0m[1mn[22m[0m[1mk[22mDeficientException low[0m[1mr[22m[0m[1ma[22m[0m[1mn[22m[0m[1mk[22mupdate low[0m[1mr[22m[0m[1ma[22m[0m[1mn[22m[0m[1mk[22mupdate! low[0m[1mr[22m[0m[1ma[22m[0m[1mn[22m[0m[1mk[22mdowndate
    





```
rank(A::AbstractMatrix; atol::Real=0, rtol::Real=atol>0 ? 0 : n*ϵ)
rank(A::AbstractMatrix, rtol::Real)
```

Compute the rank of a matrix by counting how many singular values of `A` have magnitude greater than `max(atol, rtol*σ₁)` where `σ₁` is `A`'s largest singular value. `atol` and `rtol` are the absolute and relative tolerances, respectively. The default relative tolerance is `n*ϵ`, where `n` is the size of the smallest dimension of `A`, and `ϵ` is the [`eps`](@ref) of the element type of `A`.

!!! compat "Julia 1.1"
    The `atol` and `rtol` keyword arguments requires at least Julia 1.1. In Julia 1.0 `rtol` is available as a positional argument, but this will be deprecated in Julia 2.0.


# Examples

```jldoctest
julia> rank(Matrix(I, 3, 3))
3

julia> rank(diagm(0 => [1, 0, 2]))
2

julia> rank(diagm(0 => [1, 0.001, 2]), rtol=0.1)
2

julia> rank(diagm(0 => [1, 0.001, 2]), rtol=0.00001)
3

julia> rank(diagm(0 => [1, 0.001, 2]), atol=1.5)
1
```




## 3. RAPPELS IMPORTANTS

Soit le problème linéaire suivant :
\begin{align*}
\text{Min } c^t x\\
\text{Sujet à } Ax = b\\
x ≥ 0
\end{align*}

$A$ une matrice de $m$ lignes et $n$ colonnes $(m \times n)$ et $b$ un vecteur colonne de $m$ lignes.
On suppose $m\leq n$ et $A$ est supposée de rang plein (i.e. les lignes de $A$ sont linéairement indépendantes, $rg(A) = m$, ou de manière équivalente, il y a $m$ colonnes linéairement indépendantes).

$\textbf{Définition 1}$
Une sous matrice $B$ de $A$ est dite $\textbf{base}$ si $B$ est une sous-matrice carrée (c'est-à-dire formée de $m$ colonnes de $A$) inversible (ie $B^{-1}$ existe), de dimensions $m \times m$.

$\textbf{Définition 2}$
Les **variables de base** $x_B$ sont les variables associées aux colonnes de la base $B$.

$\textbf{Définition 3}$
Une solution de base associée à la base $B$, notée $w$ correspond à poser les variables hors bases à zéro (elles sont au nombre de $n-m$), et à déterminer le vecteur des variables de base $x_B = B^{-1}b$.
Une solution de base est réalisable si $x_B \geq 0$.

### Exemple 3 page 25
\begin{align}
    \text{Min} & -2x_1 - x_2 \\
    \text{Sujet à } & x_1 + \frac{8}{3}x_2 + x_3 = 4 \\
            & x_1+x_2 + x_4 = 2 \\
            & 2x_1 + x_5 = 3 \\
            & x_1, x_2, x_3, x_4, x_5 \geq 0
\end{align}

On écrit
$$
A = \begin{pmatrix}
      1 & \frac{8}{3} & 1 & 0 & 0 \\
      1 & 1           & 0 & 1 & 0 \\
      2 & 0           & 0& 0 & 1
\end{pmatrix}
$$

Enumérez à l'aide de Julia les différentes bases et déduisez la solution optimale de ce problème.


```julia
A = [1 8/3 1 0 0; 1 1 0 1 0; 2 0 0 0 1]
b = [4; 2; 3]
c = [ -2; -1 ]
```




    2-element Vector{Int64}:
     -2
     -1




```julia
dims = size(A)
m = dims[1]
n = dims[2]
```




    5



Pour connaître le nombre possible de choix de $m$ colonnes parmi $n$, nous utilisons


```julia
binomial(n,m)
```




    10



Énumérons ces 10 combinaisons possibles:
$$
A_1 = \begin{pmatrix}
      1 & \frac{8}{3} & 1  \\
      1 & 1           & 0  \\
      2 & 0           & 0
\end{pmatrix}
\quad
A_2 = \begin{pmatrix}
      1 & \frac{8}{3} & 0  \\
      1 & 1           & 1  \\
      2 & 0           & 0
\end{pmatrix}
\quad
A_3 = \begin{pmatrix}
      1 & \frac{8}{3} & 0 \\
      1 & 1           & 0 \\
      2 & 0           & 1
\end{pmatrix}
\quad
A_4 = \begin{pmatrix}
      1 &  1 & 0  \\
      1 &  0 & 1 \\
      2 &  0 & 0
\end{pmatrix}
$$
$$
A_5 = \begin{pmatrix}
      1 & 1 & 0 \\
      1 & 0 & 0 \\
      2 & 0 & 1
\end{pmatrix}
\quad
A_6 = \begin{pmatrix}
      1 & 0 & 0 \\
      1 & 1 & 0 \\
      2 & 0 & 1
\end{pmatrix}
\quad
A_7 = \begin{pmatrix}
      \frac{8}{3} &  1 & 0 \\
      1           &  0 & 1  \\
      0           &  0 & 0
\end{pmatrix}
\quad
A_8 = \begin{pmatrix}
       \frac{8}{3} & 1 & 0  \\
       1           & 0 & 0  \\
       0           & 0 & 1
\end{pmatrix}
$$
$$
A_9 = \begin{pmatrix}
       \frac{8}{3} & 0 & 0 \\
       1           & 1 & 0 \\
       0           & 0 & 1
\end{pmatrix}
\quad
A_{10} = \begin{pmatrix}
      1 & 0 & 0 \\
      0 & 1 & 0 \\
      0 & 0 & 1
\end{pmatrix}
$$



Une première approche, naïve, consiste à écrire


```julia
A1 = [ A[:,1] A[:,2] A[:,3] ];
A2 = [ A[:,1] A[:,2] A[:,4] ];
A3 = [ A[:,1] A[:,2] A[:,5] ];
A4 = [ A[:,1] A[:,3] A[:,4] ];
A5 = [ A[:,1] A[:,3] A[:,5] ];
A6 = [ A[:,1] A[:,4] A[:,5] ];
A7 = [ A[:,2] A[:,3] A[:,4] ];
A8 = [ A[:,2] A[:,3] A[:,5] ];
A9 = [ A[:,2] A[:,4] A[:,5] ];
A10 = [ A[:,3] A[:,4] A[:,5] ];

Ai=[A1, A2, A3, A4, A5, A6, A7, A8, A9, A10];
```

et puis pour chaque base :

  1. on calcule le rang de chaque base

  
  2. on résout **si possible** $y = B^{-1}b$ où $B$ est l'une des bases $Ai$ décrites plus haut $i\in \{ 1,..., 10 \} $.


```julia
for i=1:binomial(n,m)
    if rank(Ai[i]) == m
        y = Ai[i] \ b;
        println(i, ". ", Ai[i], " ", y)
    end
end
```

    1. [1.0 2.6666666666666665 1.0; 1.0 1.0 0.0; 2.0 0.0 0.0] [1.5, 0.5, 1.1666666666666667]
    2. [1.0 2.6666666666666665 0.0; 1.0 1.0 1.0; 2.0 0.0 0.0] [1.5, 0.9375, -0.4375]
    3. [1.0 2.6666666666666665 0.0; 1.0 1.0 0.0; 2.0 0.0 1.0] [0.8, 1.2000000000000002, 1.4]
    4. [1.0 1.0 0.0; 1.0 0.0 1.0; 2.0 0.0 0.0] [1.5, 2.5, 0.5]
    5. [1.0 1.0 0.0; 1.0 0.0 0.0; 2.0 0.0 1.0] [2.0, 2.0, -1.0]
    6. [1.0 0.0 0.0; 1.0 1.0 0.0; 2.0 0.0 1.0] [4.0, -2.0, -5.0]
    8. [2.6666666666666665 1.0 0.0; 1.0 0.0 0.0; 0.0 0.0 1.0] [2.0, -1.3333333333333333, 3.0]
    9. [2.6666666666666665 0.0 0.0; 1.0 1.0 0.0; 0.0 0.0 1.0] [1.5, 0.5, 3.0]
    10. [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0] [4.0, 2.0, 3.0]


Les bases réalisables sont : A1, A3, A4, A9 et A10.
Les solutions réalisables associées à ces bases sont respectivement

\begin{align*}
\begin{pmatrix} 1.5 & 0.5 & 0 & 0 & \frac{7}{6} \end{pmatrix} \\
\begin{pmatrix} 0.8 & 1.2 & 0 & 0 & 1.4 \end{pmatrix} \\
\begin{pmatrix} 1.5 & 0 & 2.5 & 0 & 0.5 \end{pmatrix} \\
\begin{pmatrix} 0 & 1.5 & 0 & 0.5 & 3 \end{pmatrix} \\
\begin{pmatrix} 0 & 0 & 4 & 2 & 3 \end{pmatrix} \\
\end{align*}

Finalement, on compare les valeurs de fonction économique pour trouver la solution optimale.


```julia
costs = []
x = [ 1.5 ; 0.5 ]
costs = vcat( costs, dot(c, x) )
x = [ 0.8 ; 1.2 ]
costs = vcat( costs, dot(c, x) )
x = [ 1.5 ; 0 ]
costs = vcat( costs, dot(c, x) )
x = [ 0 ; 1.5 ]
costs = vcat( costs, dot(c, x) ) # c' * x
x = [ 0 ; 0 ]
costs = vcat( costs, dot(c, x) )
costs
```




    5-element Vector{Any}:
     -3.5
     -2.8
     -3.0
     -1.5
      0



La valeur optimale vaut donc -3.5, avec comme solution optimale associée $x_1 = 1.5$, $x_2 = 0.5$.

Évidemment, cette manière de procéder manque d'élégance.

Julia dispose de nombreuses librairies utiles. Nous pouvons ainsi identifier les combinaisons de colonnes avec la commande *combinations* de la librairie *Combinatorics*.

Pour ajouter ce package, il faut faire
```
using Pkg
Pkg.add("Combinatorics")
```


```julia
using Combinatorics
```

Nous pouvons calculer toutes les combinaisons possibles de colonnes avec la commande


```julia
comb = collect(combinations(1:n,m))
```




    10-element Vector{Vector{Int64}}:
     [1, 2, 3]
     [1, 2, 4]
     [1, 2, 5]
     [1, 3, 4]
     [1, 3, 5]
     [1, 4, 5]
     [2, 3, 4]
     [2, 3, 5]
     [2, 4, 5]
     [3, 4, 5]



Nous obtenons alors la $i^e$ sous-matrice avec la commande


```julia
i = 1
A[:,comb[i]]
```




    3×3 Matrix{Float64}:
     1.0  2.66667  1.0
     1.0  1.0      0.0
     2.0  0.0      0.0



Ainsi, pour énumérer les bases, nous pouvons écrire


```julia
bases = []
for i = 1:length(comb)
    B = A[:,comb[i]]
    if (rank(B) == m)
        bases = vcat(bases, i)
    end
end
bases
```




    9-element Vector{Any}:
      1
      2
      3
      4
      5
      6
      8
      9
     10



Limitons-nous aux bases réalisables


```julia
bases = []
sols = []
for i = 1:length(comb)
    B = A[:,comb[i]]
    if (rank(B) == m)
        y = B\b # équivaut mathématiquement à y = B^{-1}*b
        if all( y .>= 0)
            bases = vcat(bases, i)
            sol = zeros(n)
            sol[comb[i]] = y
            sols = vcat(sols, [sol])
        end
    end
end
sols
```




    5-element Vector{Any}:
     [1.5, 0.5, 1.1666666666666667, 0.0, 0.0]
     [0.8, 1.2000000000000002, 0.0, 0.0, 1.4]
     [1.5, 0.0, 2.5, 0.5, 0.0]
     [0.0, 1.5, 0.0, 0.5, 3.0]
     [0.0, 0.0, 4.0, 2.0, 3.0]



Nous pouvons ainsi chercher la solution optimale en cherchant celle qui donne la plus petite valeur de la fonction objectif.


```julia
idx = 1
opt = c'*sols[1][1:2]
for i = 2:length(sols)
    temp = c'*sols[i][1:2]
    if (temp < opt)
        idx = i
        opt = temp
    end
end
```

La solution et la valeur optimales sont donc respectivement


```julia
sols[idx], opt
```




    ([1.5, 0.5, 1.1666666666666667, 0.0, 0.0], -3.5)



## 4. Application

Source: Marcel Bodgan (2018), "Multiple solutions in linear programming problem", Procedia Manufacturing 22, pp 1063-1068

Nous considérons un problème de fourniture électrique. Supposons qu'il y a deux consommateurs $C_1$, $C_2$ et deux sources $R_1$, $R_2$. Les puissances nécessaires à couvrir pour les consommateurs sont $P_1$ et $P_2$, respectivement. Les distances sont $l_{ij}$, $i \in \lbrace 1,2\rbrace$, $j \in \lbrace 1,2 \rbrace$, de la source $S_i$ au consommateur $C_j$. Les puissances inconnues à transporter sont $x1=P11$, $x2=P21$, $x3=P12$, $x4=P22$. La fonction à optimiser est
$$
l11x1+l21x2+l12x3+l22x4.
$$
Avec les valeurs nominales des données, la forme canonique est
3x1+3x2+2x3+2x4−→min−x1−x2≤−10−x3−x4≤−20x1,x2,x3,x4≥


```julia

```

9614397b682d6865a1fd02e2c7717392fb2b2c35
ipynb
Jupyter Notebook
LP/.ipynb_checkpoints/Demo1-checkpoint.ipynb
fbastin/optim
07a202fc3af1e061a5580a0308949699b15a39a7
MIT
LP/.ipynb_checkpoints/Demo1-checkpoint.ipynb
fbastin/optim
07a202fc3af1e061a5580a0308949699b15a39a7
MIT
LP/.ipynb_checkpoints/Demo1-checkpoint.ipynb
fbastin/optim
07a202fc3af1e061a5580a0308949699b15a39a7
MIT
Qwen/Qwen-72B
1. YES
2. YES