[Automazione] codice matlab predator-prey
Sto provando ad implementare un algoritmo matlab per sviluppare un sistemo non lineare a tempo discreto ,che controlli la popolazione predatori e prede.
il mio sistema di partenza è:
Nella Task1 della traccia mi chiede di : "Implement the discrete-time dynamics in Matlab. Explore the behavior of the system by giving constant inputs."
Io ho iniziato scrivendo il codice in questo modo , ma non sono sicura e volevo avere un suggerimento.. e poi volevo chiedere se il sistema devo linealizzarlo da subito o basta faro dopo ,perchè nella task2 successiva mi chiede questo : "Suppose to have constant set-point populations xp,⋆, xs,⋆ (can be freely chosen). Design
an infinite horizon linear quadratic regulator to track those populations (with cost weights can be freely chosen). Notice that you need to linearize the system dynamics about the reference set-point. Simulate the evolution of the controlled system for different initial conditions."
il mio sistema di partenza è:
Nella Task1 della traccia mi chiede di : "Implement the discrete-time dynamics in Matlab. Explore the behavior of the system by giving constant inputs."
Io ho iniziato scrivendo il codice in questo modo , ma non sono sicura e volevo avere un suggerimento.. e poi volevo chiedere se il sistema devo linealizzarlo da subito o basta faro dopo ,perchè nella task2 successiva mi chiede questo : "Suppose to have constant set-point populations xp,⋆, xs,⋆ (can be freely chosen). Design
an infinite horizon linear quadratic regulator to track those populations (with cost weights can be freely chosen). Notice that you need to linearize the system dynamics about the reference set-point. Simulate the evolution of the controlled system for different initial conditions."
function TASK1 %parametri N1= 3;%identifica la dimensione del vettore colonna di ingresso ukp=rand(N1,1)*4;%vettore colonna di crescita prede con intervallo 0,4 uks= rand(N1,1)*1;%vettore colonna di decrescita predatori con intervallo 0,1 rs= 5e-6;%fattore crescita predatori numero reale rp= 1;%fattore decrescita prede numero reale %numero di stati Kmax=3; %iniziallizzare xkp=rand(1,N1)*1;%vettori di stati delle prede e poi predatori in un intervallo tra 0,1 xks=rand(1,N1)*1; %xkp(1)=0; %xks(1)=0; %formula for k=1:(Kmax-1) xkp(k+1)=(ukp(k)-uks(k))*xkp(k)*(1-xkp(k)); xks(k+1)=rs*xks(k)*(1-xks(k))*(1-uks(k)); end %plot Xvals=1:Kmax; plot(Xvals,xkp,'b',Xvals,xks,'r') end
Risposte
Alcuni suggerimenti
1) Imporre che comunque il numero delle prede e il numero dei predatori calcolati ad ogni passo siano numeri non negativi
2) I valori degli ingressi, costanti durante tutto il run, va bene che siano scelti randomicamente all'inizio, ma i valori stessi dovrebbero essere disponibili (scritti o plottati). Lo scopo infatti è capire come gli ingressi influenzano l'evoluzione del sistema
3) stesso come sopra per le condizioni iniziali.
4) Aumentare il valore di Kmax. Se si riesce a raggiungere un regime, tipicamente oscillatorio, questo avviene dopo parecchi step.
La linearizzazione serve nella task2. Le matrici che A e B che descrivono il sistema linearizzato sono utilizzate nella definizione del controllore ottimo.
Spero di esserti stato utile.
1) Imporre che comunque il numero delle prede e il numero dei predatori calcolati ad ogni passo siano numeri non negativi
2) I valori degli ingressi, costanti durante tutto il run, va bene che siano scelti randomicamente all'inizio, ma i valori stessi dovrebbero essere disponibili (scritti o plottati). Lo scopo infatti è capire come gli ingressi influenzano l'evoluzione del sistema
3) stesso come sopra per le condizioni iniziali.
4) Aumentare il valore di Kmax. Se si riesce a raggiungere un regime, tipicamente oscillatorio, questo avviene dopo parecchi step.
La linearizzazione serve nella task2. Le matrici che A e B che descrivono il sistema linearizzato sono utilizzate nella definizione del controllore ottimo.
Spero di esserti stato utile.
Grazie mille, molto utili i suggerimenti. posso chiederti gentilmente se mi dici come sviluppare il secondo punto che mi hai detto.
Dal codice mi pare che up e us siano anch'essi dei vettori per cui dovresti riuscire ad es. a fare il plot come per xp e xs.
Il plot saranno al momento 2 linee piatte ma comunque nella task 2 il trend sarà necessario perchè in tal caso, da quello che capisco, up e us dovrebbero essere le azioni di comando in uscita dal controllore ottimo, azioni che non saranno costanti ma che varieranno con k.
Nella task1 lo scopo, come detto, è quello di capire l'effetto sul processo (azione in anello aperto) al variare di tali parametri e questo permetterà di verificare se il controllore ottimo, una volta collegato al processo, si comporta in modo sensato (oppure se c'è qualcosa di sbagliato).
La cosa migliore sarebbe, durante i run per variazioni di up e us , mantenere le stesse condizioni iniziali per tutti i run, perchè se continua a variare sempre tutto non si capisce nulla, e quindi ti suggerirei di fare le prove senza randomizzare le condizioni iniziali, ma per alcune condizioni iniziali fissate di test, cosa che viene molto bene anche per le verifiche di task2.
Il plot saranno al momento 2 linee piatte ma comunque nella task 2 il trend sarà necessario perchè in tal caso, da quello che capisco, up e us dovrebbero essere le azioni di comando in uscita dal controllore ottimo, azioni che non saranno costanti ma che varieranno con k.
Nella task1 lo scopo, come detto, è quello di capire l'effetto sul processo (azione in anello aperto) al variare di tali parametri e questo permetterà di verificare se il controllore ottimo, una volta collegato al processo, si comporta in modo sensato (oppure se c'è qualcosa di sbagliato).
La cosa migliore sarebbe, durante i run per variazioni di up e us , mantenere le stesse condizioni iniziali per tutti i run, perchè se continua a variare sempre tutto non si capisce nulla, e quindi ti suggerirei di fare le prove senza randomizzare le condizioni iniziali, ma per alcune condizioni iniziali fissate di test, cosa che viene molto bene anche per le verifiche di task2.
Grazie... un'ultima domanda: sto ridefinendo i vettori up e us e volevo chiedere se ,secondo voi, il fatto che c'è scritto che up appartiene a valori tra 0 e 4 ,significa che gli elementi dei vettori possono assumere solo valori reali tra zero e 4.
Ed eventualmente, come si scriverebbe un vettore del genere in Matlab?
Ed eventualmente, come si scriverebbe un vettore del genere in Matlab?
Come li hai ridefiniti esattamente e/o come vorresti che vengano up e us?
L'ho definito così:
xx=[0.6,1,0.5,0.2;0.6,1,0.3,0.4];
uu=[2,2,0.5,0.5,2,1,1,4,1,1,1,1; 0.01,0.01,0.02,0.03,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05];
rs= 1;
for k =1:maxK-1 % p=1 s=2
xx(1,k+1)=(uu(1,k)-xx(2,k))*xx(1,k)*(1-xx(1,k));
xx(2,k+1)=rs*xx(2,k)*(1-xx(2,k))*xx(1,k)*(1-uu(2,k));
end
fig=figure;
stairs((1:maxK),[xx(1,1:maxK);xx(2,1:maxK)]','o-','LineWidth',1);
Nella seconda task , mi chiede di progettare un controllore lineare quadratico a orizzonte infinito... Come faccio a scrivere su matlab una formula con un orizzonte infinito? Non va in loop il programma?
e poi... come si fa a passare da un sistema non lineare a uno lineare su matlab?
e poi... come si fa a passare da un sistema non lineare a uno lineare su matlab?
Il codice lo devo controllare, ma per quello che riguarda il controllore:
1) Linearizzazione
Il testo ti dice innanzitutto di scegliere dei valori di set point costanti per xp e xs.
Poi, nel caso continuo mi sembra che la funzione “trim” calcoli il punto di equilibrio e la funzione “linmod” calcoli le matrici del sistema linearizzato attorno al punto di equilibrio specificato (X,U). Non so se sono utilizzabili per il caso discreto o ne esistano di equivalenti.
Per i punto di equilibrio, si ammette che X ovvero xp e xs siano pari al valore di set point. Per la U ovvero per up e us prova a mettere qualche valore desunto dalle prove precedenti (la soluzione banale con up=uk=0 non è possibile se non con xp=xs=0) e poi varia fino a che ottieni un regime stabile per i valori di X scelti.
2) Controllore Ottimo
Il fatto che il controllo sia su orizzonte infinito permette di ottenere il controllore come soluzione dell'equazione algebrica di Riccati e non come soluzione dell'equazione differenziale di Riccati. Quindi non devi fare nessun loop infinito.
Devi invece preliminarmente aver linearizzato e definito dei pesi Q e R rispettivamente per gli stati X (ovvero delle penalità per gli scostamenti eccessivi dal setpoint) e per le azioni di comando U (ovvero delle penalità per le azioni di comando eccessive). Dopodiché puoi impostare l'equazione di calcolo del regolatore (probabilmente già disponibile nel toolbox dedicato al controllo avanzato) e che comunque puoi trovare qui sia per il caso continuo che discontinuo.
http://www.dii.unimore.it/~lbiagiotti/M ... Ottimo.pdf
1) Linearizzazione
Il testo ti dice innanzitutto di scegliere dei valori di set point costanti per xp e xs.
Poi, nel caso continuo mi sembra che la funzione “trim” calcoli il punto di equilibrio e la funzione “linmod” calcoli le matrici del sistema linearizzato attorno al punto di equilibrio specificato (X,U). Non so se sono utilizzabili per il caso discreto o ne esistano di equivalenti.
Per i punto di equilibrio, si ammette che X ovvero xp e xs siano pari al valore di set point. Per la U ovvero per up e us prova a mettere qualche valore desunto dalle prove precedenti (la soluzione banale con up=uk=0 non è possibile se non con xp=xs=0) e poi varia fino a che ottieni un regime stabile per i valori di X scelti.
2) Controllore Ottimo
Il fatto che il controllo sia su orizzonte infinito permette di ottenere il controllore come soluzione dell'equazione algebrica di Riccati e non come soluzione dell'equazione differenziale di Riccati. Quindi non devi fare nessun loop infinito.
Devi invece preliminarmente aver linearizzato e definito dei pesi Q e R rispettivamente per gli stati X (ovvero delle penalità per gli scostamenti eccessivi dal setpoint) e per le azioni di comando U (ovvero delle penalità per le azioni di comando eccessive). Dopodiché puoi impostare l'equazione di calcolo del regolatore (probabilmente già disponibile nel toolbox dedicato al controllo avanzato) e che comunque puoi trovare qui sia per il caso continuo che discontinuo.
http://www.dii.unimore.it/~lbiagiotti/M ... Ottimo.pdf
Sto provando a definire i punti d'equilibrio per passare da un sistema non lineare a uno lineare. Da quello che ho capito dovrebbero uscire due coppie d'equilibrio, una per i predatori e una per le prede.
$(x_(e)^p, u_(e)^p) $ e $(x_(e)^s, u_(e)^s) $
Per trovarle sto imponendo che le due equazioni siano uguali a zero:
Nel primo caso:
$ (u_(k)^p - x_(k)^s)(x_(k)^p - (x_(k)^p)^2)=0 $ e ottengo che questa è uguale a zero se e solo se :
$ u_(k)^p= x_(k)^s $ e $ x_(k)^p= 1 $ e $ x_(k)^p= 0 $
Nel secondo caso :
$ r^s *x_(k)^s (1 - x_(k)^s)x_(k)^p (1- u_(k)^s)=0$ e ottengo che questa è uguale a zero se e solo se :
$ u_(k)^s= 1 $ e $ x_(k)^p= 0 $ e $ x_(k)^s= 1 $ e $ x_(k)^s= 0 $
Se fosse possibile vorrei avere un riscontro da parte vostra perchè ,essendo che devo avere coppie di valori, non ho capito come scegliere tra i valori diversi ottenuti dalla stessa variabile...
E' corretto come l'ho svolto per arrivare ai punti di equilibrio?
$(x_(e)^p, u_(e)^p) $ e $(x_(e)^s, u_(e)^s) $
Per trovarle sto imponendo che le due equazioni siano uguali a zero:
Nel primo caso:
$ (u_(k)^p - x_(k)^s)(x_(k)^p - (x_(k)^p)^2)=0 $ e ottengo che questa è uguale a zero se e solo se :
$ u_(k)^p= x_(k)^s $ e $ x_(k)^p= 1 $ e $ x_(k)^p= 0 $
Nel secondo caso :
$ r^s *x_(k)^s (1 - x_(k)^s)x_(k)^p (1- u_(k)^s)=0$ e ottengo che questa è uguale a zero se e solo se :
$ u_(k)^s= 1 $ e $ x_(k)^p= 0 $ e $ x_(k)^s= 1 $ e $ x_(k)^s= 0 $
Se fosse possibile vorrei avere un riscontro da parte vostra perchè ,essendo che devo avere coppie di valori, non ho capito come scegliere tra i valori diversi ottenuti dalla stessa variabile...
E' corretto come l'ho svolto per arrivare ai punti di equilibrio?
Questo è un sistema discreto e quindi in condizioni di regime stazionario deve risultare che le equazioni non devono essere zero ma fornire un punto fisso, ovvero che il valore k+1 sia uguale al valore k. Quindi:
$x^p =(u^p-x^s)*x^p*(1-x^p)$
$x^s = r^s*x^s*(1-x^s)*x^p*(1-u^s)$
Nota che ho tolto il pedice "k" per quanto detto sopra.
La soluzione $x^p=0$ e $x^s =0$ non si considera accettabile. Quindi deve valere:
$(u^p-x^s)*(1-x^p)=1$
$r^s*(1-x^s)*x^p*(1-u^s)=1$
Assegnando a $x^p$ e $x^s$ i valori di setpoint scelti, si possono trovare le u corrispondenti (attenzione che le soluzioni stiano nel campo di variazione fornito per le u e, nel caso, aggiustati il parametro r perchè lo siano) e quindi si può linearizzare.
$x^p =(u^p-x^s)*x^p*(1-x^p)$
$x^s = r^s*x^s*(1-x^s)*x^p*(1-u^s)$
Nota che ho tolto il pedice "k" per quanto detto sopra.
La soluzione $x^p=0$ e $x^s =0$ non si considera accettabile. Quindi deve valere:
$(u^p-x^s)*(1-x^p)=1$
$r^s*(1-x^s)*x^p*(1-u^s)=1$
Assegnando a $x^p$ e $x^s$ i valori di setpoint scelti, si possono trovare le u corrispondenti (attenzione che le soluzioni stiano nel campo di variazione fornito per le u e, nel caso, aggiustati il parametro r perchè lo siano) e quindi si può linearizzare.
"ingres":
$(u^p-x^s)*(1-x^p)=1$
$r^s*(1-x^s)*x^p*(1-u^s)=1$
Ho capito che in condizioni stazionarie non deve essere uguale a zero, ma perchè uguale a 1? posso mettere qualsiasi valore al di sopra di 0 ?
Non ho posto uguale ad 1 ma semplificato l'espressione precedente dividendo per $x^p$ i due membri della prima equazione e per $x^s$ i due membri della seconda equazione.
La semplificazione è lecita perchè non le soluzioni nulle non sono di interesse.
La semplificazione è lecita perchè non le soluzioni nulle non sono di interesse.
Ma io non ho capito che intende per valori di set point. Perché devo inserire questi valori per trovare la coppia di equilibrio?
A parte che non è possibile trovare 4 valori da un sistema di 2 equazioni, stiamo parlando di un sistema dinamico controllato ovvero di un sistema che deve raggiungere dei valori finali (set-point)che decidi tu.
Quindi, in base al testo "Suppose to have constant set-point populations xp,⋆, xs,⋆ (can be freely chosen)", scegli due valori costanti xp* e xs*. Il controllore che andrai a progettare dovrà comparare questi valori con i valori di xp e xs e generare dei valori up e us tali da condurre il sistema ad annullare la differenza e a regime avere xp=xp* e xs=xs*.
Pertanto possiamo ipotizzare che il sistema lavori prossimo ai valori di regime xp* e xs* e alle rispettive azioni di comando up* e us*, che per quanto detto:
1) rappresentano un punto di regime costante ovvero di equilibrio
2) rappresentano la scelta naturale del punto intorno al quale linearizzare
La scelta di xp* e xs* è a tuo piacimento, purchè all'interno del range ammesso. Risolvendo il sistema troverai up* e us* e potrai usare il parametro r per evitare di uscire dal range di ammissibilità.
Esempio: scelti xp*=0.6, xs*=0.2 , dalla prima equazione si trova up*=2.7 e sostituendo nella seconda si trova us*=1-2.0833/r e quindi per stare dentro il range di us scelgo r=4, ottenendo us*=0.4792.
Quindi, in base al testo "Suppose to have constant set-point populations xp,⋆, xs,⋆ (can be freely chosen)", scegli due valori costanti xp* e xs*. Il controllore che andrai a progettare dovrà comparare questi valori con i valori di xp e xs e generare dei valori up e us tali da condurre il sistema ad annullare la differenza e a regime avere xp=xp* e xs=xs*.
Pertanto possiamo ipotizzare che il sistema lavori prossimo ai valori di regime xp* e xs* e alle rispettive azioni di comando up* e us*, che per quanto detto:
1) rappresentano un punto di regime costante ovvero di equilibrio
2) rappresentano la scelta naturale del punto intorno al quale linearizzare
La scelta di xp* e xs* è a tuo piacimento, purchè all'interno del range ammesso. Risolvendo il sistema troverai up* e us* e potrai usare il parametro r per evitare di uscire dal range di ammissibilità.
Esempio: scelti xp*=0.6, xs*=0.2 , dalla prima equazione si trova up*=2.7 e sostituendo nella seconda si trova us*=1-2.0833/r e quindi per stare dentro il range di us scelgo r=4, ottenendo us*=0.4792.
Grazie , mi stai aiutando molto...
una volta trovati i valori della matrice A e B per linearizzare il sistema, ho riscritto il codice nel seguente modo e ho inserito anche la parte di LQR. Quando avvio il programma, mi dice che c'è un errore nel plot ma non riesco a capire perchè..
Mi esce questa scritta:
una volta trovati i valori della matrice A e B per linearizzare il sistema, ho riscritto il codice nel seguente modo e ho inserito anche la parte di LQR. Quando avvio il programma, mi dice che c'è un errore nel plot ma non riesco a capire perchè..
close all; clc; clear all;
N=5; maxK=10;
x=[0.6,1,0.5,0.2;
0.6,1,0.3,0.4];
u=[0.5,2,0.5,0.5,2,1,1,4,1,2,;
0.2,0.01,0.02,0.03,0.5,0.3,0.05,0.4,0.05,0.05];
rs= 1;
for k =1:maxK-1 % p=1 s=2
x(1,k+1)=(u(1,k)-x(2,k))*x(1,k)*(1-x(1,k));
x(2,k+1)=rs*x(2,k)*(1-x(2,k))*x(1,k)*(1-u(2,k));
end
fig=figure;
stairs((1:maxK),[x(1,1:maxK);x(2,1:maxK)]','o-','LineWidth',1);
%task2
A=[-0.5, 0.456;
0.667,0.333];
B=[0.144 , 0.6;
0.4 , 0];
Q = eye(2);
QN = Q;
rho = 1;
R = rho;
N = 5;
P = zeros(2,2,N);
K = zeros(2,2,N);% matrice di zeri 2xN
P(:,:,end) = QN;
for k = N-1:-1:1
Pkp = P(:,:,k+1);
Pk = Q + A'*inv(eye(2) + Pkp*B*inv(R)*B)*Pkp*A;
P(:,:,k) = Pk;
end
for k = 1:N-1
Pkp = P(:,:,k+1);
Kk = -inv(R+B*Pkp*B)*B*Pkp*A;
uk = Kk*x(:,k);
K(:,:,k) = Kk;
u(:,k) = uk;
x(:,k+1) = A*x(:,k)+B*uk;
end
figure
plot((1:N),x,'b-')
figure
plot((1:N),u,'r-')Mi esce questa scritta:
Vectors must be the same length.
Error in prova_progetto_2 (line 51)
plot((1:N),x,'b-')
Provo a capirci e ti dico.
Grazie, mi saresti molto d'aiuto.
x, u sono matrici 2x10, mentre N=5.
Se metti 2*N scompare l'errore e ti plotta 2 tracce (prima traccia = prima riga, seconda traccia = seconda riga).
Se metti 2*N scompare l'errore e ti plotta 2 tracce (prima traccia = prima riga, seconda traccia = seconda riga).
Grazie del consiglio, ha funzionato ma il problema era anche che avevo calcolato male i valori delle matrici A e B, che mi servivano per linearizzare il sistema.
Sistemata la task 1 e la task2 , ora sto riscontrando problemi con la task3, che mi chiede :
Supponiamo di tenere traccia delle popolazioni di riferimento variabili nel tempo $x_k ^(p\ast)$ e $x_k^(s\ast)$
per tutti i $k$. Scegli l'opzione di riferimento. Progettare ad ogni passo di tempo k un controllore MPC a orizzonte finito per risolvere il problema di controllo. È necessario linearizzare la dinamica del sistema in base al riferimento variabile nel tempo. La lunghezza dell'orizzonte temporale e i pesi possono essere scelti liberamente.
Simulare l'evoluzione del sistema controllato per diverse condizioni iniziali.
Non riesco a capire se per risolvere questa task 3 devo ricalcolarmi le matrici A e B per linearizzare il sistema o se, essendo che già nella task 2 ho preso i valori di set-point per trovarmi le matrici , posso usare quelle stesse matrici anche per la task 2
Sistemata la task 1 e la task2 , ora sto riscontrando problemi con la task3, che mi chiede :
Supponiamo di tenere traccia delle popolazioni di riferimento variabili nel tempo $x_k ^(p\ast)$ e $x_k^(s\ast)$
per tutti i $k$. Scegli l'opzione di riferimento. Progettare ad ogni passo di tempo k un controllore MPC a orizzonte finito per risolvere il problema di controllo. È necessario linearizzare la dinamica del sistema in base al riferimento variabile nel tempo. La lunghezza dell'orizzonte temporale e i pesi possono essere scelti liberamente.
Simulare l'evoluzione del sistema controllato per diverse condizioni iniziali.
Non riesco a capire se per risolvere questa task 3 devo ricalcolarmi le matrici A e B per linearizzare il sistema o se, essendo che già nella task 2 ho preso i valori di set-point per trovarmi le matrici , posso usare quelle stesse matrici anche per la task 2
Ritengo che chieda di:
0) scegliere delle leggi di variazione (lenta) per xp*, xs*
1) scegliere delle condizioni iniziali
2) fornire il valore di xp*, xs* al passo k-simo
3) calcolare us*, up* di conseguenza. Per semplicità puoi mantenere costante r (vuol dire che ci saranno dei limiti nei valori di xp* e xs* che possono essere forniti ovvero nelle leggi di cui al punto 0)
4) calcolare le nuove matrici linearizzate per i valori di xp*, xs* scelti al passo k-simo e i valori calcolati di us*, up*
5) sulla base delle nuove matrici linearizzate, ricalcolare i parametri del controllore e inserirli in linea al posto di quelli precedenti
6) memorizzare l'evoluzione del sistema.
Ritornare a 1)
0) scegliere delle leggi di variazione (lenta) per xp*, xs*
1) scegliere delle condizioni iniziali
2) fornire il valore di xp*, xs* al passo k-simo
3) calcolare us*, up* di conseguenza. Per semplicità puoi mantenere costante r (vuol dire che ci saranno dei limiti nei valori di xp* e xs* che possono essere forniti ovvero nelle leggi di cui al punto 0)
4) calcolare le nuove matrici linearizzate per i valori di xp*, xs* scelti al passo k-simo e i valori calcolati di us*, up*
5) sulla base delle nuove matrici linearizzate, ricalcolare i parametri del controllore e inserirli in linea al posto di quelli precedenti
6) memorizzare l'evoluzione del sistema.
Ritornare a 1)
Accedi a tutti gli appunti
Tutor AI: studia meglio e in meno tempo