Corso di Teoria dei Sistemi
A.A. 2002/2003
Esercitazione Matlab n.1
Studio di un sistema massa – molla
Equazioni fisiche che descrivono il sistema
h
M
k
F
Molla………………………… Fm = − k ⋅ x
Pistone (smorzatore)………… Fp = − h ⋅ x&
Bilancio di forze…….. M ⋅ &x& = F + Fm + F p
x
Scelta delle variabili di stato e di ingresso
 x1

 x2
u

=
=
Riferimenti istruzioni e funzioni Matlab
definizione di variabili e matrici (e vettori)
M = 1;
k = 0.25;
h = 1;
A = [0 1; -k/M –h/M];
x
x&
= F
Equazioni differenziali
 x&1

 x&
 2
definizione di un modello state-space
sys = ss(A,B,C,D);
= x2
autovalori di una matrice
eig(A)
1
k
h
= − ⋅ x1 − ⋅ x2 + ⋅ u
M
M
M
Equazione di uscita
definizione di un vettore equispaziato
t = linspace(0,0.5,20);
t = 0:0.025:0.5;
y = x1
visualizzazione grafica del movimento libero
initial(sys,x)
Si vuole rappresentare il sistema come
visualizzazione grafica del movimento forzato
risposta ad impulso di ampiezza unitaria
impulse(sys)
 x& = A ⋅ x + B ⋅ u

y = C ⋅ x + D ⋅u
risposta a scalino di ampiezza unitaria
step(sys)
Nel caso in esame
salvataggio dei dati in vettori
(specificando il vettore tempo)
[yl,tl,xl] = initial(sys,x,t);
[yf,tf,xf] = impulse(sys,t);
 d  x1 
1   x1   0 
 0

=

k
h  ⋅   +  1  ⋅u
 dt  x2 
− M − M   x2   M 

 x1 
y
[
]
=
1
0
⋅
 x  + [0]⋅ u

 2

tracciamento di grafici (con griglia)
plot(t,yl+yf), grid on
ciclo for
for i=1:10
istruzioni del ciclo
end
Caso non lineare
 x&1


 x& 2
= x2
Utilizzo della DEE in Simulink
dee
(da riga di comando Matlab)
k ⋅ e − x1
h
1
= − 0
⋅ x1 −
⋅ x2 +
⋅u
M
M
M
1
Codice Matlab
% TEORIA DEI SISTEMI N.O.
% Esercitazione Matlab n. 1 – 11/13.03.2003
% ======================================================================
%
%
%
%
Matlab è un interprete di comandi. Gli stessi comandi che si possono
digitare dal prompt possono essere editati in un file con estensione
.m: tutte le istruzioni vengono interpretate ed eseguite in sequenza
digitando il nome del file dal prompt (senza estensione .m).
% N.B. NON ASSEGNARE MAI AD UN FILE O AD UNA VARIABILE UN NOME CHE
% INIZI CON UN CARATTERE NUMERICO!!!
% Il simbolo "%" indica i commenti
% L'istruzione per pulire il workspace cancellando tutte le variabili
% in memoria è la seguente
clear all
% SISTEMA DINAMICO RAPPRESENTATO COME "STATE-SPACE"
% ======================================================================
% Inseriamo (definiamo) le variabili (il ";" a fine riga impedisce la
% visualizzazione).
% In questo caso si tratta dei parametri del sistema in esame
M = 1.0;
k = 0.25;
h = 0.5;
% kg
% N/m
% N*s/m
% Inseriamo poi le matrici del sistema dinamico in esame: ogni
% sequenza di valori separata da un ";" rappresenta una riga della
% matrice (ad es. B è un vettore colonna, mentre C è un vettore riga)
A
B
C
D
=
=
=
=
[0 1; -k/M -h/M];
[0; 1/M];
[1 0];
[0];
%
%
%
%
%
%
%
Definiamo poi il sistema dinamico studiato come una "variabile" qui
denominata "sys". "ss" è una funzione Matlab ("state space").
Per una completa descrizione dell'utilizzo di ciascuna funzione
si digiti "help nomedellafunzione". "help" da solo restituisce un
elenco di tutti i "toolbox" (insiemi di istruzioni e funzioni)
disponibili -> "help nometoolbox" offre un elenco delle istruzioni
e funzioni del singolo toolbox
sys = ss(A,B,C,D);
% Si provi ora a lanciare "sys" dal prompt.
% MOVIMENTO LIBERO
% ======================================================================
% Visualizziamo graficamente l'andamento dell'uscita del sistema
% conseguente ad una perturbazione dello stato ("movimento libero").
% Il sistema massa - molla è all'equilibrio in x = [0 0], u = 0.
initial(sys,[1 0])
% Come specificato dall'help ("help initial"), è possibile far
2
%
%
%
%
%
%
%
%
%
%
%
calcolare alla funzione i valori dell'andamento di stato e uscita,
anziché tracciarli graficamente. Questo si può fare successivamente
con la funzione "plot".
Salviamo gli andamenti di stato e uscita, più il vettore tempo,
in Y,T,X (la sequenza delle uscite dipende da come è stata scritta
la funzione initial).
In questo esempio si pone x1 = 1, x2 = 0 (è come se il carrello
venisse spostato in x1 = 1, e poi rilasciato con velocità iniziale
nulla).
[Y,T,X] = initial(sys,[1 0]);
% si può anche specificare l'istante di tempo finale (tf)
tf = 40;
[Y,T,X] = initial(sys,[1 0],tf);
% oppure tutto il vettore tempo: equispaziato fra 0 e 40 secondi,
% e di 1000 elementi
t = linspace(0,40,1000);
[Y,T,X] = initial(sys,[1 0],t);
% (in questo caso T coinciderà con t);
% Tracciamo ora i grafici:
% 1) andamento dell'uscita (Y)
plot(T,Y), grid on
% (ascisse = tempo; ordinate = spostamento)
% ("grid on" è un'opzione che fa tracciare una griglia sul grafico)
% 2) andamento dei due stati
%
(X è una matrice con due vettori colonna corrispondenti agli
%
andamenti dei due stati x1 e x2)
plot(T,X), grid on
% (ascisse = tempo; ordinate = spostamento e velocità)
% - movimento dello stato nello spazio delle variabili
plot(X(:,1),X(:,2)), grid on
% (ascisse e ordinate = variabili di stato x1 e x2)
plot(X(:,1),X(:,2)), grid on, axis([-1 1 –1 1])
% (“axis” consente di fissare la scala della figura).
% N.B. Sia V una matrice di tot righe e colonne.
% - V(i,j) seleziona l'elemento sulla riga i e colonna j di V;
% - V(:,j) seleziona tutta la colonna j-esima;
% - V(i,:) seleziona tutta la riga i-esima;
% - V(i,end) seleziona l'ultimo elemento della riga i;
% - V(end,j) seleziona l'ultimo elemento della colonna j;
% - V(:,end) seleziona tutta l'ultima colonna;
% - V(end,:) seleziona tutta l'ultima riga.
%
%
%
%
Si provi ora a ripetere la stessa prova ponendo nei dati iniziali
del sistema h=0 (modificando il file .m e rilanciandolo; N.B. le
modifiche risultano attive dopo che il file è stato salvato).
Si noti che le oscillazioni scompaiono per valori di h >= 1.
% Per aprire più figure
figure(1)
plot(T,X), grid on
figure(2)
3
plot(X(:,1),X(:,2)), grid on
% CICLO FOR
% ======================================================================
% Ripetiamo ora quanto visto in precedenza in un ciclo for in cui
% facciamo variare il valore di h, assegnandogli di volta in volta
% un valore preso da un vettore H che definiamo
H = [0 0.25 0.5 0.75 1];
% Ad ogni esecuzione del ciclo, cambiando h andrà ridefinita la
% matrice A ed il sistema sys
for i = 1:5
h = H(i);
A = [0 1; -k/M -h/M];
sys = ss(A,B,C,D);
[Y,T,X] = initial(sys,[1 0],t);
figure(1)
plot(T,Y), grid on, hold on
figure(2)
plot(X(:,1),X(:,2)), grid on, hold on
pause
end
%
%
%
%
- il comando "hold on" consente di sovrascrivere un grafico senza
cancellare il contenuto già presente;
- il comando "pause" attende che l'utente dia "invio" per proseguire.
(N.B. "invio" va dato con Matlab finestra attiva)
% RISPOSTA AD IMPULSO (MOVIMENTO FORZATO)
% ======================================================================
% Ripristino del sistema iniziale
h = 0.5;
A = [0 1; -k/M -h/M];
sys = ss(A,B,C,D);
% Risposta ad impulso di ampiezza unitaria
impulse(sys)
% oppure
[Y,T,X] = impulse(sys,t);
%
%
%
%
%
Come nel caso precedente, tracciare gli andamenti di Y e di X.
Si noti che la risposta ad impulso tracciata coincide con quella
che si ottiene studiando il movimento libero dello stato con
condizione iniziale [0 1].
% TRENO DI IMPULSI
% ======================================================================
%
%
%
%
%
%
%
%
%
Supponiamo di sollecitare il sistema con un treno di impulsi dati
con frequenza preimpostata (ad. es uno al secondo). Inizialmente il
sistema sarà all'equilibrio, con x = [0 0]. Se dopo un secondo dal
primo impulso ricevuto lo stato si troverà in un punto differente,
conseguentemente al secondo impulso la sua evoluzione nel tempo sarà
data dalla somma di due contributi (essendo il sistema lineare, vale
il principio di sovrapposizione degli effetti):
- il movimento libero (che si calcola con la initial);
- il movimento forzato (che si calcola con la impulse);
4
%
% Si tratterà di calcolare di volta in volta gli andamenti di stato e
% uscita del sistema, concatenando i risultati in vettori Y, X e T.
t = linspace(0,1,20);
x = [0 0];
Y = [];
X = [];
T = [];
for i = 1:15
[yl,tl,xl] = initial(sys,x,t);
[yf,tf,xf] = impulse(sys,t);
Y = [Y; yl+yf];
X = [X; xl+xf];
T = [T t+(i-1)];
end
x = xf(end,:) + xl(end,:);
% N.B. t è un vettore riga, yl, yf, xl ed xf sono vettori colonna.
%
% DOMANDA: nel codice c'è un elemento ridondante che potrebbe essere
% portato fuori dal ciclo for: quale?
figure(1)
plot(T,Y), grid on
figure(2)
plot(T,X), grid on
figure(3)
plot(X(:,1),X(:,2)), grid on
% "close all" chiude tutte le figure
5