Ecco uno script matlab per modellare un sistema PID, ovvero con un controllo Proporzionale-Derivativo-Integratore, nel quale è presente un disturbo.

Viene utilizzato l’integratore ode45 sulla dinamica complessiva del sistema esteso.

File pid_function.m (contiene la funzione chiamata dalla ode)

function zdot = pid_function(t,z)

global A B K D d Q R % importa le var. globali

%K = [kp kd ki];
[K,P] = lqr(A,B,Q,R); % ricava K

%K(1,3)= 0; % Annulla il termine integrale

u = -K*z;             % ricava u da usare riga sotto
%udot=theta
%zdot= [xdot; udot]

%Dinamica complessiva del sistema esteso
zdot = A*z + B*u + D*d;

end

Il file pid_ode.m invece è il file principale da richiamare per eseguire lo scirpt.

clc           %pulisce il workspace
clear all     %canc variabili
%close all     %chiude finestre e refresha
% commentato per permettere l’hold on delle figure

global A B K D d Q R
%——————————-
% Costanti
I = 1;
A = [0 1 0; 0 0 0 ; 1 0 0 ];
B = [0;1/I; 0];
D = [0; 1; 0];

q1 = 1;
q2 = 1;
q3 = 1;
r = 1;
Q = [q1 0 0;0 q2 0; 0 0 q3];
R = [r];

%——————————
% Intervallo di integrazione
tin = 0;
tfin = 15;
tspan = [tin tfin];

%——————————
% Valori iniziali

d = 10; % disturbo costante
x0 = [2 1];
u0 = [0];
z0=[x0′;u0′];

%——————————

esegui_ode=1;  % variabile booleana di debug per ode su zdot
if (esegui_ode==1)

options = odeset(‘abstol’,1e-10,’reltol’,1e-5); %tolleranzaassoluta=10^-10 reltol=10^-5;
[t,z]=ode45(‘pid_function’,tspan,z0,options);

%Calcolo del vettore controllo u
for i=1:length(t)
% u=-K*x
u(i,:)=-K*(z(i,1:3))’;
end

end

disegna_figure=1; % variabile booleana di debug per la stampa delle figure
if (disegna_figure==1)

figure(1);
set(1, ‘NextPlot’, ‘add’)          % imposta i parametri dell’immagine
set(1,’Units’, ‘normalized’)
set(1,’Position’, [.2 0 0.6 0.8])
set(1,’NumberTitle’, ‘off’)
set(1,’Name’, ‘Spostamento e velocità angolare, controllo ed errori’)

subplot(3,1,1); %Spost. Angolare
plot(t,z(:,1))
hold all
grid on
xlabel(‘Tempo’,’fontsize’,8,’fontname’,’tahoma’,’fontweight’,’bold’)
ylabel(‘Spost. Angolare’,’fontsize’,8,’fontname’,’tahoma’,’fontweight’,’bold’)

subplot(3,1,2); %Velocità Angolare
plot(t,z(:,2))
hold all
grid on
xlabel(‘Tempo’,’fontsize’,8,’fontname’,’tahoma’,’fontweight’,’bold’)
ylabel(‘Vel, Angolare’,’fontsize’,8,’fontname’,’tahoma’,’fontweight’,’bold’)

subplot(3,1,3); %Vettore Controllo
plot(t,u)
hold all
grid on
xlabel(‘Tempo’,’fontsize’,8,’fontname’,’tahoma’,’fontweight’,’bold’)
ylabel(‘Controllo’,’fontsize’,8,’fontname’,’tahoma’,’fontweight’,’bold’)

end

K  % => K =    2.4142    2.4142    1.0000

disp(‘Script eseguito.’)

OSSERVAZIONI

Come si vede dal grafico seguente, lo spostamento angolare nel sistema P.I.D. è
nullo (curve blu e azzurra), mentre nel sistema PD, senza integratore (curve
rossa e verde) si ha uno sfasamento residuo theta_finale= 4.142 = D/Kp (il sistema pd è ottenuto impostando a zero la K dell’integratore nel file pid_function.m, decommentando la riga 8).
Cambiando le condizioni iniziali sulla x0 i risultati non cambiano.

Grafico blu :      PID con x0=[0,0]
Grafico azzurro :  PID con x0=[2,1]

Grafico verde :    PD  con x0=[2,1]
Grafico rosso :    PD  con x0=[0,0]

K =    2.4142    2.4142    1.0000