[Generico] Consigli di sistemi non lineari da studiare e controllare

[MrChopin]

Salve ragazzi devo fare un progetto di dimanica e controllo non lineare. Devo controllare un sistema non lineare con vari tipologia di controllo e avevo intenzione di fare un modello con il DeepMPC. Avete da consigliarmi qualche modello non lineare da controllare? E avete da passarmi qualche materiale di DeepMPC? Preferibilemte qualche modello non lineare con al massimo 3 variabili di stato? Anche script di matlab toolbox da consigliare ecc...?
Non posso portare questi modelli perchè bannati:
• brusselator,
• chemostato,
• oscillatore di Van der Pol,
• modello di FitzHugh-Nagumo,
• modello preda-predatore (Lotka-Volterra),
• levitatore magnetico.
• pendolo planare

[ingres]

In rete puoi trovare modelli di reali apparati da controllare di vario tipo e natura. Di seguito ne riporto alcuni che, magari ancora un pò semplificati, possono essere un buono spunto.

1) Modello di un reattore di fermentazione
https://www.aidic.it/cet/16/50/029.pdf

2) Modello di un CSTR
https://en.wikipedia.org/wiki/Continuou ... nk_reactor

3) Modello temperatura vapore surriscaldato di un power plant
https://www.researchgate.net/publicatio ... ower_Plant

4) Turbina Eolica
https://www.osti.gov/etdeweb/servlets/purl/20675096

[MrChopin]

"ingres":In rete puoi trovare modelli di reali apparati da controllare di vario tipo e natura. Di seguito ne riporto alcuni che, magari ancora un pò semplificati, possono essere un buono spunto.
....
2) Modello di un CSTR
https://en.wikipedia.org/wiki/Continuou ... nk_reactor

Eh lo so! Stavo provando a fare da due giorni il CSTR, perchè c'erano tool e materiale di matlab su DeepMPC di cui ero molto interessato soprattutto per il campo BCI attivo, e utilizzando il tool e il materiale di matalab avrei modo di confrontare i miei risultati con quelli matlab per capire se sono giusti o meno.
Ma sto impazzendo perchè non riesco a trovare i punti di equilibrio perchè le equazioni del sistema sono un casino. Sai darmi una mano per qualche comando che posso usare ho provato di tutto.



Queste sono le equazioni del sistema.
$ { ( (dC_A(t))/dt=F/V(C_(Af)(t)-C_A(t))-k_0e^(-E/(RT(t)))C_A(t) ),( (dT(t))/dt=F/V(T_f(t)-T(t))-(DeltaH)/(rhoC_p)k_0e^(-E/(RT(t)))C_A(t)-(UA)/(rhoC_p)(T(t)-T_c(t)) ):} $

Soggetto ai seguenti vincoli:

$ { ( 0<=C_A(t)<=10),( 298<=T(t)<=800 ):} $

Dove:
1) $C_A$: è la concentrazione del reagente $A$ misurato in $(kmol)/m^3$ ho considerato che la disponibilità di reagente per questa reazione chimita sia compreso tra le $0$ e le $10(kmol)/m^3$
2) $T$: è la temperatura del reattore CSTR che ho considerato compresa tra i $298 K$ altrimenti in natura non avvengono che io sappia reazioni chimiche e $800K$ (a $1000k$ è la temperatura del plasma quindi penso che questo sia un range fattibile)
3) $C_(Af)$: è la concentrazione del reagente $A$ che inserisco dall'esterno, misurato in $(kmol)/m^3$ ho considerato che la disponibilità di reagente per questa reazione chimica sia compreso tra le $0$ e le $10$ $(kmol)/m^3$ che sarebbe legato all'ingresso $u_1$ della figura che potrei fissare ad un massimo di $10(kmol)/m^3$
4) $T_(f)$: è la temperature del bocchettone dove inserisco il fluido che definisce la temperatura del reagente che sarebbe legato all'ingresso $u_2$ che posso fissare a $298K$
5) $T_(c)$: è la temperature del bocchettone dove inserisco il fluido che definisce la temperatura del reagente che sarebbe legato all'ingresso $u_3$ che posso far variare per raffredare per scambio termico il reattore CSTR

Mi sono basato su queste due fonti:
https://it.mathworks.com/help/mpc/gs/cstr-model.html
https://www.rpi.edu/dept/chem-eng/WWW/faculty/bequette/education/links_mods/cstr/

Tutto bello tutto giusto ma non riesco a fare l'analisi delle fasi ne tanto meno riesco a trovare la soluzione di queste equazioni per trovare i punti di equilibrio:

$ { ( F/V(C_(Af)(t)-C_A(t))-k_0e^(-E/(RT(t)))C_A(t)=0),( F/V(T_f(t)-T(t))-(DeltaH)/(rhoC_p)k_0e^(-E/(RT(t)))C_A(t)-(UA)/(rhoC_p)(T(t)-T_c(t))=0 ):} $

Soggetto ai seguenti vincoli:

$ { ( 0<=C_A(t)<=10),( 298<=T(t)<=800 ):} $

ho usato questi comandi Matlab:

F = 1; % Tasso di flusso volumetrico m^3/h
V = 1; % Volume reattore m^3 
R = 1.985875; % Costante dei gas ideali di Boltzmann kcal/(kmol·K)
dH = 5960; % Entalpia kcal/kmol	
E = 11843; % Energia di Reazione kcal/kmol	
k= 34930800; % Costante non termica 1/h
rho_cp = 500; % Densita moltiplicata per la capacità termica
UA = 150; % Scambio termico alla superficie
T_f = 300; % Temperatura del flusso di alimentazione in ingresso K
C_af=10; % Concentrazione del reagente A nel flusso di alimentazione in ingresso, misurata in kmol/m^3
T_c = 292;

a=F/V;
b=E/R;
c=dH*k/(rho_cp);
d=UA/(V*rho_cp);
syms x_1 x_2
equ1= (a)*(C_af-x_1)-x_1*k*exp(-(b/x_2)) ==0;
equ2= (a)*(T_f-x_2)-(c)*x_1*k*exp(-(b/x_2))-(d)*(x_2) ==0;
equ3= x_1 >=0;
equ4= x_2 >=298;
equ5= x_1 <=10;
equ6= x_2 <=800;

set_equ=[equ1, equ2, equ4, equ5, equ6];
[x_1 x_2]=vpasolve(set_equ,[x_1,x_2])

OR

[x_1 x_2]=solve(set_equ,[x_1,x_2])

Come la metto metto il risultato è lo stesso:



Niente. E per non parlare di pplane8 che va a casaccio che una volta mi caccia un punto di equilibrio altre volte no capisco le variazione parametriche possono far cambiare la topologia dei punti di equilibrio ma non capisco che errori faccio.

Una volta mi esce questo:


Una volta mi esce questo

[MrChopin]

"ingres":In rete puoi trovare modelli di reali apparati da controllare di vario tipo e natura. Di seguito ne riporto alcuni che, magari ancora un pò semplificati, possono essere un buono spunto.

1) Modello di un reattore di fermentazione
https://www.aidic.it/cet/16/50/029.pdf

2) Modello di un CSTR
https://en.wikipedia.org/wiki/Continuou ... nk_reactor

3) Modello temperatura vapore surriscaldato di un power plant
https://www.researchgate.net/publicatio ... ower_Plant

4) Turbina Eolica
https://www.osti.gov/etdeweb/servlets/purl/20675096

Ma penso per i motivi di cui sopra mi arrendo al fatto che non ce la farò mai a risolvere quelle equazioni perchè non credo che siano risolvibili e quindi che ne pensi di questo modello?

https://it.mathworks.com/help/deeplearn ... ulink.html

[ingres]

Ti dovessi dire mi sembra che manchi un T_c nella seconda equazione.

Comunque miei suggerimenti

1) Usa tranquillamente un modello più semplice. Il secondo modello di CSTR che hai individuato va altrettanto bene per lo scopo che ti proponi.

2) Il problema di non riuscire a trovare un punto di equilibrio dipende da tanti fattori. Scartate le situazioni banali di errori nel modello, può essere un problema che si parte troppo da lontano oppure delle limitazioni eccessive sugli ingressi oppure un'instabilità propria del modello, ecc.
Se i tool falliscono e si vuole proprio addivenire ad una soluzione, solitamente conviene impostare anche il problema del controllo in modo semplificato, individuando le variabili da controlllare e mettendo qualche semplice PID come controllore per generare le azioni di comando ovvero gli ingressi esterni al modello, lasciando agli stessi, almeno inizialmente, un range ampio di variazione.
Si decidono dei set realistici e guadagni cautelativi per i PID e quindi si attiva la simulazione e si vede se si riesce a convergere ad una situazione di regime, che è il punto di equilibrio cercato.

Fammi sapere se sei riuscito ad andare avanti

[MrChopin]

"ingres":Ti dovessi dire mi sembra che manchi un T_c nella seconda equazione.

Comunque miei suggerimenti

1) Usa tranquillamente un modello più semplice. Il secondo modello di CSTR che hai individuato va altrettanto bene per lo scopo che ti proponi.

2) Il problema di non riuscire a trovare un punto di equilibrio dipende da tanti fattori. Scartate le situazioni banali di errori nel modello, può essere un problema che si parte troppo da lontano oppure delle limitazioni eccessive sugli ingressi oppure un'instabilità propria del modello, ecc.
Se i tool falliscono e si vuole proprio addivenire ad una soluzione, solitamente conviene impostare anche il problema del controllo in modo semplificato, individuando le variabili da controlllare e mettendo qualche semplice PID come controllore per generare le azioni di comando ovvero gli ingressi esterni al modello, lasciando agli stessi, almeno inizialmente, un range ampio di variazione.
Si decidono dei set realistici e guadagni cautelativi per i PID e quindi si attiva la simulazione e si vede se si riesce a convergere ad una situazione di regime, che è il punto di equilibrio cercato.

Fammi sapere se sei riuscito ad andare avanti

Ti faccio fare due risate ho chiesto a chatgpt usando una vpn mi ha stampato questo codice:

% Definizione delle variabili simboliche
syms x_1 x_2

% Definizione dei parametri noti
a = 1; % valore noto di a
C_af = 0.5; % valore noto di C_af
T_f = 400; % valore noto di T_f
k = 0.1; % valore noto di k
b = 0.2; % valore noto di b
c = 0.3; % valore noto di c
d = 0.4; % valore noto di d

% Definizione delle equazioni
equ1 = a*(C_af - x_1) - x_1*k*exp(-b/x_2);
equ2 = a*(T_f - x_2) - c*x_1*k*exp(-b/x_2) - d*x_2;
equ3 = x_1 >= 0;
equ4 = x_2 >= 298;
equ5 = x_1 <= 10;
equ6 = x_2 <= 800;

% Definizione della funzione di obiettivo
fun = @(x) double([isAlways(equ1); isAlways(equ2); isAlways(equ3); isAlways(equ4); isAlways(equ5); isAlways(equ6)]);

% Risoluzione del sistema di equazioni
x0 = [0; 400]; % punto iniziale della soluzione
sol = fsolve(fun, x0);

% Estrazione dei risultati
x_1_sol = sol(1);
x_2_sol = sol(2);

% Stampa dei risultati
fprintf('Risultati:\n');
fprintf('x_1: %.4f\n', x_1_sol);
fprintf('x_2: %.4f\n', x_2_sol);

Usa fsolve per risolverla ma come risultato mi esce che il punto di equilibrio stia a x1=0 ovvero assenza del concetrato A e x2 tempeuratura a 400K gli ho fatto fare anche il codice matlab per il DeepMPC (Spolier non funziona)

[ingres]

Che sono gli stessi valori che mi sembra assuma come punto iniziale di innesco per il calcolo.
Alla fine non è stato in grado di trovare la soluzione e ha mantenuto i valori iniziali.

Comunque CA = 9.5 e T=293.57 risulterebbero essere i valori di equilibrio. Prova a ricontrollare le equazioni e a rilassare il vincolo inferiore su T.

[MrChopin]

"ingres": Comunque CA = 9.5 e T=293.57 risulterebbero essere i valori di equilibrio. Prova a ricontrollare le equazioni e a rilassare il vincolo inferiore su T.

Mi puoi passare lo script per favore al meno capisco l'errore magari mi perdo in un bicchiere d'acqua

[ingres]

Non ho in questo momento a disposizione Matlab, quindi ho banalmente riscritto le equazioni in Excel e usato il relativo Solver.

Ma se prendi le equazioni che hai già scritto in Matlab, controlli bene la seconda equazione (come ti ho detto mi sembra che manchi la sottrazione con T_c) e anzichè imporre che T>298 imponi T>290, con qualunque solver di Matlab secondo me arrivi al risultato.

Prova e sappimi dire.

[MrChopin]

"ingres":Non ho in questo momento a disposizione Matlab, quindi ho banalmente riscritto le equazioni in Excel e usato il relativo Solver.

Ma se prendi le equazioni che hai già scritto in Matlab, controlli bene la seconda equazione (come ti ho detto mi sembra che manchi la sottrazione con T_c) e anzichè imporre che T>298 imponi T>290, con qualunque solver di Matlab secondo me arrivi al risultato.

Prova e sappimi dire.

L'ho fatto ma non mi trovo con te mi trovo che

x_1 =9.999999809301705435983268588863 
x_2 = 169.6907474353380559633941211231

e purtroppo non ha senso questa soluzione perchè 170K sono -113 gradi e non ha senso una reazione chimica a quella temperature. Ma qua secondo me è una questione parametrica perchè secondo me cambia in base ai parametri. Mi passi i tuoi?

[ingres]

Gli stessi che hai inserito nel tuo esempio e che ho posto in colonna B righe da 1 a 10
F/V 1 1/h
R 1.985875 kcal/mole K
dH 5960 kcal/mole
E 11843 kcal/mole
K 34930800 1/h
rho_cp 500
UA 150
T_f 300 K
C_af 10 kmol/m3
T_c 292 K

Ho poi messo i valori di C_A e T nelle celle B12 e B13 e queste sono le equazioni:
=B1*(B9-B12)-B5*B12*EXP(-B4/B2/B13)
=B1*(B8-B13)-B3/B6*B5*B12*EXP(-B4/B2/B13)-B7/B6*(B13-B10)

[MrChopin]

la differenza sta qua secondo me B3/B6*B5*B12 fa prima la moltiplicazione e poi la divisione e quindi cambia tutto prova a cambiare questo punto e fare (B3/B6)*B5*B12 e vedi se ti trovi in linea con me

[ingres]

Fatto, ma non cambia. Excel esegue i calcoli in ordine di priorità e per la stessa priorità in sequenza.

Comunque é sospetto che ti venga una T così bassa.

Hai fatto la modifica che ti avevo indicato su
equ2 = a*(T_f - x_2) - c*x_1*k*exp(-b/x_2) - d*x_2

scrivendo la stessa così:
equ2 = a*(T_f - x_2) - c*x_1*k*exp(-b/x_2) - d*(x_2-T_c)

[MrChopin]

"ingres":Fatto, ma non cambia. Excel esegue i calcoli in ordine di priorità e per la stessa priorità in sequenza.

Comunque é sospetto che ti venga una T così bassa.

Hai fatto la modifica che ti avevo indicato su
equ2 = a*(T_f - x_2) - c*x_1*k*exp(-b/x_2) - d*x_2

scrivendo la stessa così:
equ2 = a*(T_f - x_2) - c*x_1*k*exp(-b/x_2) - d*(x_2-T_c)

Allora si trova e sto impazzendo inutilmente dietro al punto di equilibrio perchè sono un faggiano. Effettivamente a me che me frega dove sta sto punto di equilibrio? Tanto male che vada uso Tc che è la legge di contro $u_3$ e faccio variare il punto di equilibrio. Giusto? Oppure mi sbaglio?

[MrChopin]

"ingres":Gli stessi che hai inserito nel tuo esempio e che ho posto in colonna B righe da 1 a 10
F/V 1 1/h
R 1.985875 kcal/mole K
dH 5960 kcal/mole
E 11843 kcal/mole
K 34930800 1/h
rho_cp 500
UA 150
T_f 300 K
C_af 10 kmol/m3
T_c 292 K

Ho poi messo i valori di C_A e T nelle celle B12 e B13 e queste sono le equazioni:
=B1*(B9-B12)-B5*B12*EXP(-B4/B2/B13)
=B1*(B8-B13)-B3/B6*B5*B12*EXP(-B4/B2/B13)-B7/B6*(B13-B10)

Comunque avevo ragione tu fai $(DeltaH)/(rhoC_pk_0)$ infatti facendo così mi trovo i tuoi stessi valori

a=F/V;
b=E/R;
c=dH*k/(rho_cp*k);
d=UA/(V*rho_cp);

syms x_1 x_2
equ1= (a)*(C_af-x_1)-x_1*k*exp(-(b/x_2)) ==0;
equ2= (a)*(T_f-x_2)-(c)*x_1*exp(-(b/x_2))-(d)*(x_2-T_c) ==0;

Infatti poi si trova con te

[MrChopin]

"ingres":Fatto, ma non cambia. Excel esegue i calcoli in ordine di priorità e per la stessa priorità in sequenza.

Comunque é sospetto che ti venga una T così bassa.

Hai fatto la modifica che ti avevo indicato su
equ2 = a*(T_f - x_2) - c*x_1*k*exp(-b/x_2) - d*x_2

scrivendo la stessa così:
equ2 = a*(T_f - x_2) - c*x_1*k*exp(-b/x_2) - d*(x_2-T_c)

No sono il re dei faggiani ho capito l'errore in questa equazione mettevo due volte la costante integrativa

equ2= (a)*(T_f-x_2)-(c)*k*x_1*exp(-(b/x_2))-(d)*(x_2-T_c) ==0;

perchè la mettevo sia nell'equazione che nel termine c

c=dH*k/(rho_cp);

E quindi non mi trovavo grazie per lo spunto ingres sono impazzito per 5 giorni.
Mo a forza lo devo controllare sto sistema è una questione di principio!

[MrChopin]

per dovere di cronaca alla fine ce l'ho fatta ma presenta un comportamento strano cioè se io considero l'intorno del punto di equilibrio si comporta normalmente ha sto punto di equilibrio in 9.5 e 230. Ok tutto bello tutto giusto, ma se mi allontano un pò è come se avesse una sottospecie di sella in cui le curve o puntano verso il punto di equilibrio o vanno ad infinito, che però giustamente lo jacobiano non riporta. Secondo te mi posso fidare di questo piano delle fasi? Posso sfruttare anche quella sottospecie di "zona di sella" in qualche modo o molto probabilmente è solo un errore dell'ode45 dovuto ai termini di ordine superiore? Perchè più zoommo sulla "zona sella" e più si diradano le curve

[ingres]

Difficile a dirsi. I sistemi non lineari possono benissimo avere dei comportamenti strani, ma al momento non mi preoccuperei.
Il tuo scopo, se ho inteso bene, non è tanto il sistema in anello aperto, ma progettare un controllo DeepMPC. Sarà compito del sistema controllato tenere il sistema nel punto di lavoro.

[MrChopin]

Sto cercando di fare una finezza volevo cercare di risolvere l'equazione parametricamente per poter non doverla risolvere a mano l'equazione. Però mi dice che forse non è possibile ottenere una soluzione esplicita al problema come potrei fare?

Questo è lo script:

syms C_A T a b c d k C_af T_f

% equ1= diff(C_A,t)==(a)*(C_af-C_A)-C_A*k*exp(-(b/T))
% 
% equ2= diff(T,t)==(a)*(T_f-T)-(c)*C_A*exp(-(b/T))-(d)*(T)

equ1= (a)*(C_af-C_A)-C_A*k*exp(-(b/T))==0

equ2= (a)*(T_f-T)-(c)*C_A*exp(-(b/T))-(d)*(T)==0

set_equ=[equ1, equ2];

syms F V R dH E k rho C_p UA T_f C_af T_c

set_equ= subs (set_equ, [a b c d], [F/V, E/R, dH*k/(rho*C_p), UA/(V*rho*C_p)])

[solC_A,sol_T] =solve(set_equ, [C_A, T])

[ingres]

Se intendi trovare la soluzione analitica del sistema per dei generici parametri F/V, E/R, ecc. non penso che sia fattibile. Puoi trovare C_A in funzione di T e ricondurti così ad un'equazione solo in T, ma risolvere quest'ultima mi sembra impresa ardua in quanto si tratta di un'equazione trascendente.

Alternative ne vedo poche. L'unico suggerimento che mi viene in mente è che se il termine esponenziale di Arrhenius non variasse troppo nell'effettivo campo di lavoro, potresti approssimarlo con una costante (ovvero ammettiamo che la "costante di velocità" sia effettivamente costante) e in tal caso risulta un sistema lineare facilmente risolvibile.

[MrChopin]

"ingres":Se intendi trovare la soluzione analitica del sistema per dei generici parametri F/V, E/R, ecc. non penso che sia fattibile. Puoi trovare C_A in funzione di T e ricondurti così ad un'equazione solo in T, ma risolvere quest'ultima mi sembra impresa ardua in quanto si tratta di un'equazione trascendente.

Alternative ne vedo poche. L'unico suggerimento che mi viene in mente è che se il termine esponenziale di Arrhenius non variasse troppo nell'effettivo campo di lavoro, potresti approssimarlo con una costante (ovvero ammettiamo che la "costante di velocità" sia effettivamente costante) e in tal caso risulta un sistema lineare facilmente risolvibile.

Mi trovo su tutto e ci avevo pensato anche io ma poi come faccio a far una valutazione aprametrica del punto di equilibrio al variare dei parametri?