Partie 6: Contrôle de systèmes dynamiques II

Les parties sur les observateurs sont tirées du livre édité en 2008 par Denis Dochain et intitulé "Automatic Control of Bioprocesses".


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]:
Pour afficher le code python, cliquer sur le bouton:

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

Observateurs

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:

  1. l'observateur doit converger plus rapidement que la dynamique du système
  2. si $\hat{x}(0)=x(0)$ alors $\hat{x}(t)=x(t)$ pour tout $t>0$, ce qui veut dire que si on part de la bonne valeur de l'état, alors on doit pouvoir estimer parfaitement l'état aux temps suivants

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)


Saturations

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)