Réaliser une intégrale pour signal échantilloné en C

Bonjour, J’essaye de réaliser le calcul suivant: theta = (mesure_capteur - mesure_capteur_calculé) x mesure_capteur_calculé x gamma x Tech x (z/(z-1)). Il s’agit d’une commande adaptative dont l’explication est présente ici: https://fr.mathworks.com/help/physmod/sps/ref/modelreferenceadaptivecontroller.html voici le code que j’ai réalisé:

#include <stdio.h>

int main()
{
    float y = 19.9;
    float ym = 20.1;
    float gamma = 0.01;
    float Ts = 1.0;
    float integral = 0.0;
    float theta = 0.0;
    float theta_1 = 0.0;
    int i = 0;
    while(i < 1000)
    {
        theta = (y-ym)*ym*gamma;
        integral = integral +((theta+theta_1)/2.0)*1.0;
        theta_1 = theta;
        printf("integral = %f \n",integral);
        i++;
    }

    return 0;
}

voici aussi le bout de code sur un debugger en ligne: https://www.onlinegdb.com/Hy0SeKK28

Le problème que je rencontre est que le résultat de theta diverge vers des valeurs très grandes alors que normalement, avec les valeurs que j’ai choisi, la valeur de theta devrait être proche de 1. Je pense que le problème vient de mon calcul d’integral, sauriez vous comment corriger le problème ? Merci d’avance

J’aurais une question par rapport au PID : que signifie la variable z ?

Salut,

J’ai du mal à voir comment tu es arrivé à ce code-là.

Mon principal conseil est de bien reprendre tranquillement les équations dans la documentation de Matlab, voire même dans l’article cité en source si tu y as accès. À partir de ces équations (qui sont des transformées en z), on peut retrouver des équations aux différences (autrement dit des suites), qui sont traduisibles très facilement dans n’importe quel langage de programmation.

Ceci dit, j’ai quelques commentaires sur ton code.

• theta ne correspond pas à la définition de $\theta$ dans la page de Matlab, ça n’aide pas à la compréhension.
• Tu aussi une erreur de signe, le moins de la formule de la doc n’apparaît pas dans ton code.
• On dirait que tu essaies de faire la méthode du trapèze pour ton intégration. Le truc, c’est que c’est déjà un contrôleur discret, donc tu as juste à prendre les équations telles qu’elles sont écrites pour que ça marche. En voulant intégrer à ta sauce, tu risque surtout de rendre les paramètres inadaptés.
• Ton code fait (bon an, mal an) une intégration, mais c’est normal d’avoir une valeur qui augmente beaucoup, parce que tu intègres… une constante. Tu ne peux pas avoir le "vrai" comportement sans avoir la boucle de rétroaction complète.

Tu ne l’as pas expliqué : pourquoi tu veux faire ce calcul concrètement ?

Bonjour,

Salut Math’s,

La variable z correspond à la variable de la transformé en z

Salut Aabu,

Merci d’avoir pris le temps de lire en détail mon code et l’article matlab.

• Je suis désolé pour la différence de nom entre le theta de mon code et le nom utilisé dans matlab, je n’avais pas fait attention à cela.

• Merci également pour l’erreur de signe, j’ai corrigé ça dans le debugger en ligne: https://www.onlinegdb.com/r1Li1Q5nU

• Sinon j’ai effectivement effectué la methode du trapèze pour réaliser mon intégrale, je n’ai pas trouvé d’autres moyens de la réaliser sur un controleur. Je me suis également documenté plus en détail sur ce type de commande comme tu me l’avais conseillé, je n’ai pas réussi à lire la référence indiqué par matlab car il s’agit d’un livre mais j’ai trouvé ce cours qui utilise les mêmes équations et qui semble s’inspirer de la même source que matlab: http://www.phoneoximeter.org/uploads/media/EECE574–11_MRAC_01.pdf A la page 18/73 du pdf, on peut lire l’équation de θ (je suis désolé de ne pas avoir copié l’équation sur le site, je n’ai pas encore trouvé comment inclure des équations mathématiques) De ce que l’on peut y lire, j’ai l’impression que l’intégration est nécessaire.

Je n’avais pas pensé que le problème qui pourrait venir de la boucle de retroaction, je vais essayer de simuler un système plus complet sur matlab, je posterai le résultat des que je l’aurais fait.

Sinon j’essaye de réaliser ce calcul afin de pouvoir l’utiliser dans le contrôle d’un moteur à courant continu, le but de l’expérience est surtout pédagogique.

Merci à vous d’eux d’eux d’avoir pris le temps de répondre à mon message

Alors, oui l’intégration est nécessaire, parce que sinon tu ne fais plus le même calcul.

Par contre, j’ai l’impression que tu te mélanges les pinceaux entre les formalismes continus et discrets.

• La transformée en z correspond à du temps discret. Pour l’implémenter, il faut juste traduire les transformées en équations aux différences. Il n’y a rien de plus à faire et il n’y pas d’approximation d’intégrale à faire.
• La transformée de Laplace correspond à du temps continu. Pour implémenter un contrôleur continu dans un dispositif à temps discret (comme un microcontrôleur typiquement), il faut le discrétiser. Ça passe par l’approximation des intégrales, dont la méthode des trapèze fait partie, mais ce n’est pas la seule. Ça revient à implémenter une équation aux différences.

En gros le parallèle entre les deux est le suivant.

On peut en général passer de l’une à l’autre en faisant des approximations. En général, on a plus tendance à approximer le temps continu par le temps discret.

Les chemins habituels pour obtenir une équation aux différences programmable dans n’importe quel langage :

• on a déjà une transformée en z et il suffit de la transformer en équation aux différences ;
• ou alors on a une transformée de Laplace et il faut alors trouver une transformée en z approchante, puis obtenir l’équation aux différences correspondante ;
• ou alors on a une équation différentielle décrivant le système, on applique un schéma d’intégration numérique et on obtient une équation aux différences sans passer par les transformations fréquentielles (par exemple ce que je fais ici).

Merci pour les explications, j’ai essayé de mettre en pratique la seconde solution que tu as proposé, consistant à obtenir l’équation à implémenter dans le controleur en partant de la transformée de laplace. J’ai récupéré cet exemple matlab qui me parait simple à comprendre et qui est accompagné d’un pdf expliquant le fonctionnement du projet: https://fr.mathworks.com/matlabcentral/fileexchange/44416-simple-adaptive-control-example Voici une capture d’écran du fichier simulink:

Le calcul de θ me parait identique à celui que j’essaie de faire et l’intégrateur correspond à 1/p.

Le résultat de θ lorsque l’on souhaite que θ n’ait aucun effet sur la régulation est le suivant:

Pour mon controleur en temps discret j’ai donc essayé d’écrire l’intégrateur suivant:

J’ai obtenu la fonction de transfert sur scilab de 1/p pour une periode d’échantillonage de 1s grâce au bout de code suivant:

clear;

s=poly(0,'s');

sl=syslin('c',(1)/(s));   // Continuous-time system in transfer form
slss=tf2ss(sl);                  // Now in state-space form
sld=cls2dls(slss,1);          // Continuous to discrete transform
sldt=ss2tf(sld);                 // Converts in transfer form
disp(sldt);                      // Display z-transform

La fonction que j’ai obtenue a été de la forme: Y(z)/(Xz) = (0.5+0.5 * z)/(z-1) J’ai ensuite ramené cette équation sous forme d’équation de récurrence parceque je suis plus à l’aise avec cette forme que les équations aux différences que je ne connais pas:

Y(z) (z-1) = X(z) (0.5 + 0.5 * z)

Y(z) = 0.5 X(z) z^(-1) + 0.5 X(z) + Y(z) z^(-1)

Y(n) = Y(n-1) + 0.5 X(n) + 0.5 X(n-1)

Ce qui est assez marrant puisque je suis retombé sur le résultat donné par la methode du trapèze.

J’ai aussi essayé de simuler cette fonction dans scilab avec mon régulateur dans une boucle de retroaction en faisant en sorte que θ n’ait pas non plus d’effets sur la régulation mais malheureusement le résultat a une fois de plus divergé comme sur le debugger en ligne.

J’avoue que je ne comprend pas ce que j’ai fait de mal cette fois , peut être que les variables que j’ai choisi pour l’intégration n’étaient pas les bonnes ? J’ai pris Y(n) = θ, Y(n-1) = θ(n-1), X(n) = (valeur_mesuré_systeme-valeur_modèle_systeme)valeur_modèle_systeme-1.0*gamma)

Ce qui est assez marrant puisque je suis retombé sur le résultat donné par la methode du trapèze.

En fait, c’est pas très marrant, c’est parce que, d’après la doc, cls2dls fait la transformation bilinéaire, autrement dit la méthode du trapèze.

J’ai aussi essayé de simuler cette fonction dans scilab avec mon régulateur dans une boucle de retroaction en faisant en sorte que θ n’ait pas non plus d’effets sur la régulation mais malheureusement le résultat a une fois de plus divergé comme sur le debugger en ligne.

Sans le code, ça va être difficile d’aider. Ça peut dépendre vraiment de comment tu as fait tout le reste.

Rebonjour, Oui pas de soucis voici le code, j’avais juste évité de le poster car il est assez long, je pense que la partie intéressante se situe dans la boucle for:

clear;

s=poly(0,'s');

sl=syslin('c',(1)/(s));   // Continuous-time system in transfer form
slss=tf2ss(sl);                  // Now in state-space form
sld=cls2dls(slss,1);          // Continuous to discrete transform
sldt=ss2tf(sld);                 // Converts in transfer form
disp(sldt);                      // Display z-transform

// E_PFCdel
// PFC to order 1 with delay and delayed disturbance
clear; xdel(winsid());
tf=4500; w=1:1:tf; // length
tech=1; // sampling
// initialization
u=zeros(1, tf);
MV=u; sm=u; CV=u; DV=u;
//
tau=33.8; //time constant first order
am=exp(-tech/tau); bm=1-am; // model
ap=exp(-tech/tau); bp=1-ap; // process could be different
r=input("r(default:200)="); // delay for instance r = 0 or 500
if r==[] then r=200; end;
km=1; // gain model
kp=1; // process
SP=100; // setpoint
trbf=65; // tuning factor: desired closed loop time response: only tuning control
lh=1-exp(-tech*3/trbf);
cal = 0;
cal_1 = 0;
theta = 0;
theta_1 = 0;

for ii= 2+r:1:tf
    DV(ii)=0;
    // process
    if ii> 800 then DV(ii)=20; end;
    CV(ii)=CV(ii-1)*ap+bp*(MV(ii-1-r)*kp+DV(ii-r));
    // model
    // to be implemented in all PLC's
    sm(ii)=sm(ii-1)*am+bm*MV(ii-1)*km;
    // control equation
    spred=CV(ii)+sm(ii)-sm(ii-r); // delay
    cal = (CV(ii)-sm(ii-r))*sm(ii-r)*-1*0.1;
    theta = theta_1 + 0.5 * cal_1 +0.5*cal
    cal_1 = cal;
    theta_1 = theta;
    disp(theta);
    d=(SP-spred)*lh+sm(ii)*bm;
    MV(ii)=d/(km*bm);
    if MV(ii)>100 then MV(ii)=100; end // constraint on MV
   // if MV(ii)>(MV(ii-1)+2) then MV(ii)=(MV(ii-1)+2); end // constraint on MV
    //if MV(ii)<(MV(ii-1)-2) then MV(ii)=(MV(ii-1)-2); end // constraint on MV
    end

scf(0);
plot(w,MV,"b",w,CV,"r",w,DV,"m"); xlabel("sec");
a=gca(); a.grid=[1,1]; a.tight_limits="on"; a.data_bounds=[0,-0;tf,160];
title("PFC CONSTRAINT and DELAY: MV b / CV r / DV m");

Encore une fois, merci pour le temps que tu prends pour m’aider

Merci beaucoup pour l’information (désolé de vous couper)! Je vais peut être enfin réussir à comprendre ces calculs.

Salut,

Le code est difficilement compréhensible, même avec le PDF du FileExchange (je n’ai pas Matlab pour ouvrir l’autre). Les noms de variables ne correspondent pas au schéma, je n’arrive pas à savoir ce qui est vraiment le contrôleur dans tout ça.

Dans ton code, j’ai l’impression que tu n’utilises pas theta après l’avoir calculé, donc il n’aura jamais l’effet escompté. Une des deux boucles n’est pas fermée en somme.

Enfin, n’hésite pas à te poser pour réfléchir. Je te vois filer de ressource en ressource, passer d’un modèle à un autre, d’un outil à un autre… C’est pas forcément la meilleure façon d’arriver à comprendre ce que tu fais.

Salut,

Oui je suis désolé, j’ai enlevé la partie ou j’utilisais thêta (MV(ii) = theta * MV(ii) ) et j’ai oublié de la remettre avant de poster mon code. Je vais suivre ton Conseil, et essayer de me poser et de réfléchir au problème à tête reposé. Je progresserais sûrement plus vite de cette manière qu’en switchant entre plein de docs. Je posterais un nouveau message des que j’aurais du nouveau