Fondamenti di Automatica: sistemi TD e mappa logistica preda-predatore

Corso di Fondamenti di Automatica

MARCO COSTANZO
marco.costanzo@unicampania.itSistemi TD

Risposta a segnali noti: gradino e sinusoide

Analisi della risposta indiciale

• Calcoliamo la risposta indiciale
  Y(z) = G(z)
  Z
  Z-1
  bmzm + bm-1zm-1 + ... + biz + bo
  =
  z
  zn + an-1zn-1 + ... + a1z + do z - 1
• Valore iniziale e finale
  y(0) = lim Y (z) =
  Z->00
  So,
  m <n
  bm,
  m=n
• Se G(z) non ha poli in 1
  yo = lim(z-1)Y(z) = G(1)
  Z->1
• Sistemi del I ordine
  G(z) =
  u(1-p)
  z - p
  Y(z) = G(z)
  Z
  Z -1
  (z-1)(z-p)
  zu(1 - p)
  uz
  -uz
  +
  z-1 z- p
  y(k) = u(1-pk),
  k ≥ 0

1
0.8
p = 0.5
0.6
p = 0.8
y/p
p = 0.9
0.4
0.2 H
0 #
0
5
10
15
20
25
30
35
40
k

2
+
1.8
1.6
....
1.4
...
1.2
-P =- 0.9
y/p 1
- p =- 0.5
0.8
0.6
0.4
0.2
0 1
0
5
10
15
20
25
30
35
40
k
All'avvicinarsi del polo all'origine, il sistema
risponde sempre più velocemente

Analisi della risposta indiciale: Sistemi del II ordine

G(z) =
u(1-p1)(1-p2)
Z - Z1
1 - Z1
(z-P1)(z-P2)

• Nel dominio del tempo la risposta è del tipo
  y(k) = p + P1pk + P2pk,
  k ≥ 0
• L'andamento della risposta dipende dalla
  posizione dello zero oltre che dei poli
• Può esserci sovraelongazione anche con modi
  aperiodici non alternanti (poli reali in | | < 1)
• La sottoelongazione nasce quando lo zero è fuori
  dal cerchio unitario (in | | > 1)

3
2.5
2
Z1 = 0.95
Z1 = 0.9
y/μ 1.5
- Z1 = 0.5
1
Z1= 0
0.5
0
0
5
10 15
20
25
30
35
40
k

1
000
0.5
z = 0
-
z = 2
-
z1 = 1.5
-0.5
-
z1 = 1.1
O
-1
-1.5
0
10
15
20
25
30 35
40
k
Poli: p1 = 0.4 e p2 = 0.8

Analisi della risposta indiciale: Sistemi del II ordine con poli complessi

G(z) =
u(1-2p cos0+p2)z
z2- 2p cos 0z + p2

• Ha poli complessi e coniugati in
  p = petje
• Nel dominio del tempo la risposta è del tipo
  y(k) = p + 2|Q|pkcos(Ok +2Q),
  k ≥ 0
• L'andamento della risposta è di tipo pseudo-
  periodico convergente al valore di regime u
  se p < 1

1.6
1.4
1.2
1
y/p 0.8
*
0.6
0.4
0.2
0 *
0
5
10 15
20
25
30
35
40
k
Poli: p = 0.9e+jT/6

Risposta a regime sinusoidale

! Dato un sistema asintoticamente stabile con fdt G (z)

• Consideriamo un ingresso sinusoidale
  ZU sin 0
  ZU sin 0
  u(k) = U sin Ok ,
  k ≥ 0 => U(z) = 72 - 2 cos 0z + 1
  (z - eje)(z - e-je)
  Assumendo, per semplicità, poli distinti per G (z), la risposta forzata vale
  Y(z) = G(z)
  zU sin ℮
  (z - eje)(z - e-je)
  =
  i=1
  Złi
  Z - Pi
  +
  z - eje
  +
  z - e-je'
  n
  zQ
  zQ*
  con Q =
  UG(eje)
  2j
  Antitrasformando
  y(k) = Στη
  n
  ripi + 2
  i=1
  k
  UG(eje)|
  2j
  cos(0k + 2G(eje) - Tt/2) =
  i=1
  > ripk + U|G(ej0)|sin(Ok + LG(eje),
  n
  k ≥ 0
• A regime «sopravvive» solo il modo proprio dell'ingresso (sinusoidale) e quindi la risposta a regime è ancora
  sinusoidale con la stessa frequenza dell'ingresso ma ampiezza e fase modificate dalla fdt valutata in eje
  yr(k) = U|G(eje)| sin(0k + 2G(eje))

Funzione di risposta in frequenza (risposta armonica)

• È la funzione complessa di variabile reale 0 € [0, 1]
  G(eje) = G(z)
  z=ej0
• cioè la restrizione della fdt alla semicirconferenza di raggio unitario eventualmente privata dei poli di G (z) su di essa
• È una funzione periodica di periodo 21 e siccome G (z) è una funzione razionale fratta con coefficienti tutti reali, vale
  G(e-je) = G*(eje)
• Ecco perché i valori di G(eJe) per 0 € [-Tt, 0] si possono ricavare da quelli per 0 € [0, 1]
• Può essere utilizzata non solo per calcolare la risposta a segnali sinusoidali ma anche a segnali periodici sviluppabili in
  serie di Fourier e segnali non periodici trasformabili secondo Fourier (Discrete Fourier Transform - DFT)
  F[f(k)] = F(eje) =
  +00
  k =- 80
  f(k)e-jek
• Che è invertibile sotto ipotesi poco stringenti
  f(k) =
  1
  2T
  1
  -TT
  F(ej0)ejek de =
  1
  TT
  0
  IF(eje)| cos (0k + 2F (eje)) de
  Per segnali reali

Risposta a regime sinusoidale: ingresso generico

• Considerato un ingresso generico u(k) trasformabile secondo Fourier
  u(k) =
  1
  0
  -
  TT
  | |U(ej0)| cos (0k + ZU(eje)) do
• Nell'ipotesi di sistema asintoticamente stabile, per linearità, la risposta è la sovrapposizione (continua) delle risposte agli
  infiniti segnali sinusoidali
  y(k) =
  1
  π
  න
  G(eje)|U(eje)|cos(0k+ZU(ej0) + 2G(eje))de
• Ma y(k) è a sua volta trasformabile secondo Fourier e quindi
  |Y(eje)| cos (0k + LY(eje)) de
  y(k) =
• In conclusione
  Y(eje) = G(eje)U(eje)
• In uscita ci sono solo armoniche già presenti in ingresso, ciascuna scalata e sfasata per la funzione di risposta armonica
• Come per i sistemi a tempo continuo la risposta in frequenza si può rappresentare graficamente tramite
• Diagrammi di Bode, polari e di Nichols
  1
  TT
  න

Schemi a blocchi

• Le interconnessioni elementari sono le stesse di quelle dei sistemi a tempo
  continuo
• Serie
  u
  y
  u
  y
  Ga(z)
  GŁ(z)
  Ga(z)GB(z)
• Parallelo
  Ga(z)
  y
  u
  y
  Ga(z) + GB(z)
  +
  GŁ(z)
• Retroazione
  u
  +
  y
  G,(z)
  u
  Ga(z)
  y
  1 + Ga(z)GB(z)
  Gb(z)
  u
  V +

Schemi a blocchi: proprietà strutturali e stabilità

Vista l'analogia formale degli schemi a blocchi e della formula per il calcolo della fdt
data la i-s-u, valgono le stesse considerazioni sulla proprietà strutturali e la stabilità
interna e quella BIBO

Sistemi in serie

• Se uno zero di Ga (z) cancella un polo di Gb (z), il sistema serie ha una parte non raggiungibile ma
  osservabile
• Se un polo di Ga (z) è cancellato da uno zero di Gb (z), il sistema serie ha una parte non osservabile ma
  raggiungibile
• Il sistema serie è asintoticamente stabile se e solo se entrambi i sistemi connessi in serie lo sono
• La stabilità BIBO di entrambi i sistemi componenti implica la stabilità BIBO del sistema serie, ma non è
  vero il viceversa

Sistemi in parallelo

• Se Ga (z) e Gb (z) hanno poli in comune, il sistema parallelo ha una parte non raggiungibile e non
  osservabile
• Il sistema parallelo è asintoticamente stabile se e solo se entrambi i sistemi connessi in parallelo lo sono
• La stabilità BIBO di entrambi i sistemi componenti è equivalente alla stabilità BIBO del sistema parallelo

Sistemi in retroazione

• Se uno zero di Ga (z) cancella un polo di Gb (z), il sistema in retroazione ha una parte non raggiungibile e
  non osservabile
• La stabilità dei sistema in retroazione non ha una relazione diretta con la stabilità dei singoli sottosistemi
• Se ci sono cancellazioni tra zeri di Ga (z) e poli di Gb (z) non asintoticamente stabili, il sistema in
  retroazione non è asintoticamente stabile

Discretizzazione di Eulero

Oltre al modello a dati campionati, un altro modo possibile per discretizzare le
equazioni dinamiche di un modello tempo continuo è quello di sostituire la
derivate rispetto al tempo con il rapporto incrementale ottenendo
x*(k + 1) - x*(k)
Ts
~ *= Ax(t) + Bu(t)
Ottenendo:
x*(k + 1) = [I + TsA]x*(k) + TsB*u(k)
Questa volta otteniamo un'approssimazione accettabile solo per valori piccoli di
Ts

Esempi di Modelli Non Lineari

La Mappa Logistica - Modello Preda-Predatore

Popolazione di prede senza predatori

x1(k+1) = x1(k) + rx1(k) -y(x1(k))2

• X1: numero di prede ad una certa generazione
• r > 0 : tasso di riproduzione
• y > 0 : coefficiente di competizione (risorse limitate e le prede competono
  per ottenerle).
• Se y = 0 le risorse sono illimitate (si usa spesso per prede erbivore)
  Nel caso di risorse illimitate (y = 0), il modello è lineare e instabile (la
  popolazione cresce indefinitamente)

La Mappa Logistica - Modello Preda-Predatore: Popolazione di predatori senza prede

X2(k + 1) = x2(k) - ux2(k)

• X2: numero di predatori ad una certa generazione
• [ E]0,1] : tasso di diminuzione netta (mortalità meno riproduzione)
  In assenza di prede i predatori di estinguerebbero (sistema asintoticamente
  stabile)

La Mappa Logistica - Modello Preda-Predatore: Popolazione di prede e predatori

x1(k + 1) = x1(k) + rx1(k) - y(x1(k))'
2
- ax1(k)x2(k)
X2(k+1) = x2(k) - ux2(k) + bx1(k)x2(k)

• Il numero di prede diminuirà proporzionalmente al numero di incontri x1x2
  scalato di un coefficiente a
• a > 0 : Efficienza del predatore nell'uccidere una preda
• La popolazione dei predatori aumenterà proporzionalmente al numero di
  incontri scalato di un coefficiente b > 0
• b # a, numero di nuovi predatori non è uguale al numero di prede uccise

La Mappa Logistica - Modello Preda-Predatore: forma classica

Assumendo y = 0, n = 1 + r, ß = 1 - u, si ottiene la forma (classica)
x1(k+1)= (n-ax2(k))x1(k)
X2(k+1) =(B+bx1(k))x2(k)

• I punti di equilibrio sono
  Xe1 = [0] , Xe2 =
  1 - 31
  b
  n-1
  a

La Mappa Logistica - Modello Preda-Predatore: Linearizzazione

Dal momento che
df
dx
=
mn- ax2
bx2
₿ + bx1
-ax1
F
Linearizziamo intorno a xe1
8x(k + 1) =
01
8x(k)
E' instabile perché gli autovalori sono n, ß e n > 1

• Linearizziamo intorno a xe2
  8x(k + 1) =
  -r
  La
  1
  b
  -
  1
  a
  b
  8x(k)
  Che ha polinomio caratteristico z2 - 2z + 1 + ru le cui radici sono 1 ± jyru
  che hanno modulo >1 -> instabile

La Mappa Logistica - Modello Preda-Predatore: versione discretizzata

In realtà il modello presentato
x1(k + 1) = x1(k) + rx1(k) - y(x1(k))
2
- ax1(k)x2(k)
X2(k+1)= x2(k) -ux2(k) +bx1(k)x2(k)
È la versione discretizzata con il metodo di Eulero di
*1(t) =rx1(t)-y(x1(t))-ax1(t)x2(t)
*2(t) = ux2(t) +bx1(t)x2(t)
Con Ts =1
Con Ts generico sarebbe:
X1 (k + 1) = x1 (k) + (rx1 (k) - Y(x1 (k)
2
- ax1(k)x2(k) ) Ts
X2(k+1)=x2(k) +(-ux2(k)+bx1(k)x2(k))Ts

La Mappa Logistica - Modello Preda-Predatore: parametri e simulazione

*1(t) = rx1(t)-v(x1(t))2-ax1(t)x2(t)
×2(t) = 1x2(t) +bx1(t)x2(t)
Dt = 0.001; %Sampling Time per Eulero|
% parametri prede
r = 1.1; %tasso riproduzione prede
gamma = 0; % coefficiente di competizione tra le prede
a = 0.4; % efficienza del predatore nell'uccidere una preda
% parametri predatori
mu = 0.4; % tasso di diminuzione netta predatori (mortalità meno riproduzione)
b = 0.1; % nuovi predatori a seguito di uccisione preda
num_steps = round(20/Dt);
% condizioni iniziali
x1(1) = 10;
x2(1) = 1;
for k=1:num_steps
x1(k+1) =x1(k) + (r*x1(k) - gamma*(x1(k)^2) - a*x1(k)*x2(k)) * Dt;
x2(k+1) =x2(k) + (-mu*x2(k) + b*x1(k)*x2(k) ) * Dt;
t(k+1) = t(k) + Dt;
end
18
0
0
prede
predatori
16
14
12
10
8
6
4
2
0
0
5
10
15
20
25

La Mappa Logistica - Modello Preda-Predatore: parametri e simulazione 2

*1(t) = rx1(t)-v(x1(t))2-ax1(t)x2(t)
*2(t) = 1x2(t) +bx1(t)x2(t)
0
0
prede
predatori
Dt = 0.001; %Sampling Time per Eulero
% parametri prede
r = 1.1; %tasso riproduzione prede
gamma = 0; % coefficiente di competizione tra le prede
a = 0.4*10; % efficienza del predatore nell'uccidere una preda
% parametri predatori
mu = 0.4; % tasso di diminuzione netta predatori (mortalità meno riproduzione)
b = 0.1; % nuovi predatori a seguito di uccisione preda
10€
8
6
4
2
0
0
1
2
3
4
5
6