In [1]:
# -*- coding: utf-8 -*-
from IPython.display import HTML
HTML('''<script>
code_show=true;
function code_toggle() {
if (code_show){
$('div.input').hide();
} else {
$('div.input').show();
}
code_show = !code_show
}
$( document ).ready(code_toggle);
</script>
Pour afficher le code python, cliquer sur le bouton:
<button onclick="javascript:code_toggle()">Afficher code python</button>
''')
Out[1]:
In [2]:
import numpy as np
from matplotlib import pyplot as plt
from ipywidgets import interact, fixed
#from IPython.html.widgets import interact, fixed
import scipy.integrate as scint
from matplotlib import patches as pat
Les lois de commande qui sont construites (PID ou autres types) dépendent de certaines variables du système (variables d'état, d'entrée ou de sortie).
Par exemple, pour le contrôle de la concentration du susbtrat $S$ dans un réacteur continu, exemple que l'on a considéré précédemment, la commande PID:
$$ Q(t)=V\mu(S^\ast)+K_p(S^\ast-S^m(t))+K_i\int_0^t(S^\ast-S^m(s))ds+K_d \frac{d(S^\ast-S^m)}{dt}$$dépend de la concentration de substrat ou plus exactement de sa mesure $S^m$.
Or, cette quantité n'est pas forcément mesurée, et d'autant moins en temps réel. Pour appliquer la commande, il va donc falloir que l'on estime la valeur de $S$. C'est ce que vont permettre de faire les observateurs.
Le principe de l'observateur est donné sur la figure suivante:
Définition: Un observateur pour le système:
$$\left\{ \begin{array}{rcl} \dot{x}&=&f(x,u)\\ y&=&h(x,u) \end{array} \right. $$est un système d'équations de la forme:
$$ \left\{ \begin{array}{rcl} \dot{z}&=&\hat{f}(z,y,u)\\ \hat{x}&=&\hat{h}(z,y,u) \end{array} \right. $$tel que $\hat{x}$ donne une "bonne" estimation de $x$, ce que l'on peut traduire par:
$$ \lim_{t\rightarrow \infty} \left\Vert x-\hat{x} \right\Vert = 0$$Parmi tous les observateurs possibles, on va chercher ceux qui ont les deux propriétés intéressantes suivantes:
Une forme d'observateur couramment utilisée, qui peut avoir les propriétés recherchées, est la suivante:
$$ \left\{ \begin{array}{rcl} \dot{z}&=&f(z,u)+C\left[h(z,u)-y\right]\\ \hat{x}&=&z \end{array} \right. $$L'équation en $z$ est constitué d'un premier terme qui copie la dynamique du système à observer, et d'un second terme correctif qui dépend de l'erreur entre l'estimation donnée par l'observateur et la sortie (mesurée) $y$ du système.
Exemple: Croissance d'une population dans un réacteur continu (suite)
Un modèle du système est donné par:
$$ \left\{ \begin{array}{crl} \frac{dB}{dt}= & \mu(S)B &-\frac{Q}{V}B\\ \frac{dS}{dt}= & -k\mu(S)B&+\frac{Q}{V}(S_0-S) \\ \end{array} \right. $$Supposons que l'on ne mesure que le dégagement de gaz représenté par la quantité $\mu(S)B$. Alors on a: $$ y = \mu(S)B =: h(S,B)$$
Un observateur du système est donné par: $$ \left\{ \begin{array}{crl} \frac{d\hat{B}}{dt}= & \mu(\hat{S})\hat{B} &-\frac{Q}{V}\hat{B}+C_1\left[\mu(\hat{S})\hat{B}-y\right]\\ \frac{d\hat{S}}{dt}= & -k\mu(\hat{S})\hat{B}&+\frac{Q}{V}(S_0-\hat{S})+C_2\left[\mu(\hat{S})\hat{B}-y\right]\\ \end{array} \right. $$
In [3]:
def reacteur_obs(x,t,k,muast,KS,KI,Qin,V,S0,Sast,control_type,coeffcontrol,obs,Cobs,disturb):
B = x[0] #biomass
S = x[1] #substrat
if control_type in ['PI','PID']:
indObs = 3
Shat = x[indObs+1]
dx = np.zeros(6)
dx[2] = Sast-S
dx[5] = Sast-Shat
else:
dx = np.zeros(4)
indObs = 2
if obs == 0:
xest = x[0:indObs]
elif obs == 1:
xest = x[indObs:]
Q = fonction_u(t,xest,Sast,k,muast,KI,KS,V,S0,control_type,coeffcontrol)*(1+disturb)
mu = muast*S/(KS+S+S**2/KI)
dx[0] = mu*B-Q/V*B
dx[1] = -k*mu*B+Q/V*(S0-S)
Bhat = x[indObs]
Shat = x[indObs+1]
muhat = muast*Shat/(KS+Shat+Shat**2/KI)
dx[indObs] = muhat*Bhat-Q/V*Bhat+Cobs[0]*(muhat*Bhat-mu*B)
dx[indObs+1] = -k*muhat*Bhat+Q/V*(S0-Shat)+Cobs[1]*(muhat*Bhat-mu*B)
return dx
def culture_cont2(Sast,control_type,coeffcontrol,obs,coeffobs,disturb):
tmax = 100
temps = np.linspace(0,tmax,2000) # vecteur temps
k = 0.6; muast = 2.3; KS = 10; KI = 0.1; S0 = 3.2; B0 = 9; Qin = 0.01; V = 0.5;
B0obs = coeffobs[0]
S0obs = coeffobs[1]
Cobs = coeffobs[2:]
# integration numerique de l'EDO
if control_type in ['PI','PID']:
x0 = np.array([B0,S0,0,B0obs,S0obs,0])
indObs = 3
else:
x0 = np.array([B0,S0,B0obs,S0obs])
indObs = 2
x = scint.odeint(reacteur_obs,x0,temps,args=(k,muast,KS,KI,Qin,V,S0,Sast, \
control_type,coeffcontrol,obs,Cobs,disturb))
if obs == 0:
u = fonction_u(temps,x[:,0:indObs],Sast,k,muast,KI,KS,V,S0,control_type,coeffcontrol)
elif obs == 1:
u = fonction_u(temps,x[:,indObs:],Sast,k,muast,KI,KS,V,S0,control_type,coeffcontrol)
plt.figure(figsize = (10, 3))
plt.subplots_adjust(hspace=0.4,wspace=0.4)
plt.subplot2grid((1,2),(0,0))
plt.plot(temps,x[:,0],'r',label='Biomasse')
plt.plot(temps,x[:,1],'g',label='Substrat')
plt.plot(temps,x[:,indObs],'r--',label='Biomasse estimée')
plt.plot(temps,x[:,indObs+1],'g--',label='Substrat estimé')
plt.plot(np.array([0,temps[-1]]),np.array([Sast,Sast]),'k--',label='S*')
plt.legend(); plt.xlabel('time (h)')
plt.subplot2grid((1,2),(0,1))
plt.plot(temps,u,'r',label='Debit')
plt.legend(); plt.xlabel('time (h)')
plt.show()
# Contrôle boucle fermée de la culture bactérienne dans un réacteur continu: commande proportionnelle intégrale dérivée
# ---------------------------------------------------------------------------------------------------------------------
# Loi de commande
def fonction_u(t,x,Sast,k,muast,KI,KS,V,S0,control_type,coeffcontrol):
if control_type == 'BO': # BOUCLE OUVERTE
# loi donnée par Qast=mu(Sast)*V
Qast = muast*Sast/(KS+Sast+Sast**2/KI)*V
if type(t)==float: # cas où t est scalaire
valu = Qast
else: # cas où t est un vecteur
valu = np.ones(len(t))*Qast
elif control_type == 'P': # BOUCLE FERMEE action PROPORTIONNELLE
# récupération des paramètres de la loi de commande
kprop = coeffcontrol
# et de la valeur de S
if type(t)==float: # cas où t est scalaire
valS = x[1]
else: # cas où t est un vecteur
valS = x[:,1]
# loi donnée par mu(Sast)*V+kprop*(Sast-S)
valu = muast*Sast/(KS+Sast+Sast**2/KI)*V+kprop*(Sast-valS)
elif control_type == 'PI': # BOUCLE FERMEE action PROPORTIONNELLE INTEGRALE
# récupération des paramètres de la loi de commande
kprop = coeffcontrol[0]
kint = coeffcontrol[1]
# et de la valeur de valint, qui est l'intégrale entre 0 et t de Sast-S
if type(t)==float: # cas où t est scalaire
valint = x[2]; valS = x[1]
else: # cas où t est un vecteur
valint = x[:,2]; valS = x[:,1]
# loi donnée par mu(Sast)*V+kprop*(Sast-S) + kint*valint
# où valint est l'intégrale entre 0 et t de Sast-S
Qast = muast*Sast/(KS+Sast+Sast**2/KI)*V
valu = Qast+kprop*(Sast-valS)+kint*valint
elif control_type == 'PID': # BOUCLE FERMEE action PROPORTIONNELLE INTEGRALE DERIVEE
# récupération des paramètres de la loi de commande
kprop = coeffcontrol[0]
kint = coeffcontrol[1]
kderiv = coeffcontrol[2]
# et des valeurs de valint (intégrale de Sast-S), de S et de B
if type(t)==float: # cas où t est scalaire
valint = x[2]; valS = x[1]; valB=x[0]
else: # cas où t est un vecteur
valint = x[:,2]; valS = x[:,1]; valB = x[:,0]
# loi donnée par mu(Sast)*V+kprop*(Sast-S) + kint*valint + kderiv*(dSast/dt-dS/dt)
mu = muast*valS/(KS+valS+valS**2/KI)
Qast = muast*Sast/(KS+Sast+Sast**2/KI)*V
valu = (Qast+kprop*(Sast-valS)+kint*valint+kderiv*k*mu*valB)/(1+kderiv/V*(S0-valS))
return valu
Test de l'observateur, mais sans l'intégrer à la loi de commande
In [4]:
# Simulation sans utiliser les observations dans la loi de commande: on suppose l'état connu
culture_cont2(2,'PI',np.array([0.07,0.01]),0,np.array([13,4,-2,0.1]),0.2)
Test du couplage de l'observateur avec la loi de commande
In [5]:
# Simulation en utilisant les observations dans la loi de commande
culture_cont2(2,'PI',np.array([0.07,0.01]),1,np.array([13,4,-2,0.1]),0.2)
Dans le paragraphe précédent, on a vu que, lorsque certaines variables du système n'étaient pas mesurées, il était possible de les estimer grâce à un observateur.
Un autre problème que l'on peut rencontrer en pratique vient des saturations sur les commandes. En effet, pour des raisons souvent physiques, on ne peut pas appliquer n'importe quelle commande à notre système. Il y a en effet des limites physiques aux valeurs de commande, que l'on peut exprimer de la manière suivante:
$$ u_{min} ≤ u ≤ u_{max}$$Exemple: Croissance d'une population dans un réacteur continu (suite)
Dans cet exemple la commande étant un débit d'entrée, celle-ci ne pourra pas être négative. Les limites physiques de la pompe impose également une valeur de débit maximale, que l'on notera $Q_{max}$ et que l'on prendra égale à $0.1 L/h$. On a donc: $$ 0 ≤ Q ≤ Q_{max}$$ Or, on voit sur la figure précédente que la loi de commande que l'on a utilisée propose des valeurs de $Q$ négatives et supérieures à $Q_{max}$. C'est donc une loi de commande que l'on ne pourra pas appliquer telle qu'elle au système.
Un première solution à ce problème consiste à appliquer une saturation à la commande calculée $u_{calc}$.
La valeur $u_{réelle}$ que l'on appliquera sera déterminée à partir de $u_{calc}$ de la manière suivante:
$$ u_{réelle} = \left\{ \begin{array}{ll} u_{min} & \text{si }u_{calc} ≤ u_{min}\\ u_{calc} & \text{si }u_{min} < u_{calc} < u_{max}\\ u_{max} & \text{si }u_{max} ≤ u_{calc}\\ \end{array}\right.$$
In [6]:
def reacteur_obs(x,t,k,muast,KS,KI,Qin,V,S0,Qmax,Sast,control_type,coeffcontrol,obs,Cobs,sat,coeffsat,disturb):
B = x[0] #biomass
S = x[1] #substrat
if control_type in ['PI','PID']:
indObs = 3
Shat = x[indObs+1]
dx = np.zeros(6)
dx[2] = Sast-S
dx[5] = Sast-Shat
else:
dx = np.zeros(4)
indObs = 2
if obs == 0:
xest = x[0:indObs]
elif obs == 1:
xest = x[indObs:]
Q = fonction_u(t,xest,Sast,k,muast,KI,KS,V,S0,control_type,coeffcontrol)*(1+disturb)
if sat == 1:
if Q<0: Q=0
elif Q>Qmax: Q=Qmax
mu = muast*S/(KS+S+S**2/KI)
dx[0] = mu*B-Q/V*B
dx[1] = -k*mu*B+Q/V*(S0-S)
Bhat = x[indObs]
Shat = x[indObs+1]
muhat = muast*Shat/(KS+Shat+Shat**2/KI)
dx[indObs] = muhat*Bhat-Q/V*Bhat+Cobs[0]*(muhat*Bhat-mu*B)
dx[indObs+1] = -k*muhat*Bhat+Q/V*(S0-Shat)+Cobs[1]*(muhat*Bhat-mu*B)
return dx
def culture_cont3(Sast,control_type,coeffcontrol,obs,coeffobs,sat,coeffsat,disturb):
tmax = 100
temps = np.linspace(0,tmax,2000) # vecteur temps
k = 0.6; muast = 2.3; KS = 10; KI = 0.1; S0 = 3.2; B0 = 9; Qin = 0.01; V = 0.5; Qmax = 0.1;
B0obs = coeffobs[0]
S0obs = coeffobs[1]
Cobs = coeffobs[2:]
# integration numerique de l'EDO
if control_type in ['PI','PID']:
x0 = np.array([B0,S0,0,B0obs,S0obs,0])
indObs = 3
else:
x0 = np.array([B0,S0,B0obs,S0obs])
indObs = 2
x = scint.odeint(reacteur_obs,x0,temps,args=(k,muast,KS,KI,Qin,V,S0,Qmax,Sast, \
control_type,coeffcontrol,obs,Cobs,sat,coeffsat,disturb))
if obs == 0:
u = fonction_u(temps,x[:,0:indObs],Sast,k,muast,KI,KS,V,S0,control_type,coeffcontrol)
elif obs == 1:
u = fonction_u(temps,x[:,indObs:],Sast,k,muast,KI,KS,V,S0,control_type,coeffcontrol)
if sat > 0: u = (u<=0)*0+(0<u)*(u<Qmax)*u+(u>Qmax)*Qmax
plt.figure(figsize = (10, 3))
plt.subplots_adjust(hspace=0.4,wspace=0.4)
plt.subplot2grid((1,2),(0,0))
plt.plot(temps,x[:,0],'r',label='Biomasse')
plt.plot(temps,x[:,1],'g',label='Substrat')
plt.plot(temps,x[:,indObs],'r--',label='Biomasse estimée')
plt.plot(temps,x[:,indObs+1],'g--',label='Substrat estimé')
plt.plot(np.array([0,temps[-1]]),np.array([Sast,Sast]),'k--',label='S*')
plt.legend(); plt.xlabel('time (h)')
plt.subplot2grid((1,2),(0,1))
plt.plot(temps,u,'r',label='Debit')
plt.legend(); plt.xlabel('time (h)')
plt.show()
In [7]:
# Simulation avec loi de commande PI, observateur et saturation
culture_cont3(2,'PI',np.array([0.07,0.01]),1,np.array([13,4,-2,0.1]),1,0,0.2)
Imposer une saturation sur une commande peut néanmoins dégrader fortement la dynamique du système en boucle fermée, voire carrément déstabiliser le système. Cela est essentiellement dû au fait que lorsque la commande est saturée, le terme intégral continu de grossir car il continue d'intégrer l'erreur: on parle d'effet "wind-up".
Méthode anti-windup
Pour pallier ce problème, on utilise des techniques d'"anti-wind-up" qui consistent en fait à vider le terme intégral (c'est à dire le diminuer) lorsque la commande est saturée. Concrètement, au lieu de calculer le terme intégral
$$I=K_i\int_0^t(y^\ast-y^m(s))ds$$comme présenté précédemment, c'est à dire en résolvant:
$$ \frac{dI}{dt}=K_i(y^\ast-y^m),$$on va résoudre l'équation:
$$ \frac{dI}{dt}=K_i(y^\ast-y^m)+K_{aw}(u_{réel}-u_{calc})$$.
In [8]:
def reacteur_obs(x,t,k,muast,KS,KI,Qin,V,S0,Qmax,Sast,control_type,coeffcontrol,obs,Cobs,sat,coeffsat,disturb):
B = x[0] #biomass
S = x[1] #substrat
if control_type in ['PI','PID']:
indObs = 3
Shat = x[indObs+1]
dx = np.zeros(6)
dx[2] = Sast-S
dx[5] = Sast-Shat
else:
dx = np.zeros(4)
indObs = 2
if obs == 0:
xest = x[0:indObs]
elif obs == 1:
xest = x[indObs:]
Q = fonction_u(t,xest,Sast,k,muast,KI,KS,V,S0,control_type,coeffcontrol)*(1+disturb)
if sat == 1:
if Q<0: Q=0
elif Q>Qmax: Q=Qmax
elif sat == 2:
if Q<0: Qr=0
elif Q>Qmax: Qr=Qmax
else:
Qr = Q
if control_type in ['PI','PID']:
dx[2] = Sast-S+coeffsat*(Qr-Q)
dx[5] = Sast-Shat+coeffsat*(Qr-Q)
Q = Qr
mu = muast*S/(KS+S+S**2/KI)
dx[0] = mu*B-Q/V*B
dx[1] = -k*mu*B+Q/V*(S0-S)
Bhat = x[indObs]
Shat = x[indObs+1]
muhat = muast*Shat/(KS+Shat+Shat**2/KI)
dx[indObs] = muhat*Bhat-Q/V*Bhat+Cobs[0]*(muhat*Bhat-mu*B)
dx[indObs+1] = -k*muhat*Bhat+Q/V*(S0-Shat)+Cobs[1]*(muhat*Bhat-mu*B)
return dx
In [12]:
# Simulation avec loi de commande PI, observateur et saturation et anti-windup
culture_cont3(2,'PI',np.array([0.07,0.01]),1,np.array([7,4,-2,0.1]),2,100,0.2)