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

[MrChopin]

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?

[ingres]

Ci sono un paio di strade percorribili.

Se al variare dei parametri in un range assegnato il punto di equilibrio (o di funzionamento) non varia troppo, allora si possono linearizzare le equazioni attorno al punto di equilibrio nominale e quindi trovare come tale punto varia parametricamente.

Qualora ci sia una variabile principale che contraddistingue il punto di funzionamento, si può assegnare ai restanti parametri dei valori nominali, trovare la curva che descrive il punto di funzionamento al variare del parametro principale e quindi determinare le curve di correzione per tener conto degli altri parametri.

Ad esempio, se ammettiamo che il parametro principale sia T e che le altre grandezze, eccetto C_A, possano vedersi come variazioni rispetto a dei valori nominali, si può calcolare la curva C_A=C_A(T) per tali valori nominali e poi, fissandosi su un punto a T (e quindi C_A) assegnato, vedere i delta di C_A al variare degli altri parametri presi singolarmente. Si otterrano così n curve di correzione e se le variazioni sono piccole basterà sommare gli effetti di ciascuna curva.

Però hai proprio necessità di fare questa valutazione?

[MrChopin]

Ci sono un paio di strade percorribili.

Se al variare dei parametri in un range assegnato il punto di equilibrio (o di funzionamento) non varia troppo, allora si possono linearizzare le equazioni attorno al punto di equilibrio nominale e quindi trovare come tale punto varia parametricamente.

Qualora ci sia una variabile principale che contraddistingue il punto di funzionamento, si può assegnare ai restanti parametri dei valori nominali, trovare la curva che descrive il punto di funzionamento al variare del parametro principale e quindi determinare le curve di correzione per tener conto degli altri parametri.

Ad esempio, se ammettiamo che il parametro principale sia T e che le altre grandezze, eccetto C_A, possano vedersi come variazioni rispetto a dei valori nominali, si può calcolare la curva C_A=C_A(T) per tali valori nominali e poi, fissandosi su un punto a T (e quindi C_A) assegnato, vedere i delta di C_A al variare degli altri parametri presi singolarmente. Si otterrano così n curve di correzione e se le variazioni sono piccole basterà sommare gli effetti di ciascuna curva.

Però hai proprio necessità di fare questa valutazione?

Beh si le vedo anche in latri progetti per l'esame e quindi penso che sia necessario, mo ci provo e vediamo.

[MrChopin]

Comuqnue per tenerti informato ce l'ho fatta anche con lo sviluppo in serie di taylor facendo uno sviluppo in serie di taylor di entrambe le variabili intorno ai loro punti di equilibrio, vedi se ha senso logico o invece non ha senso logico. Perchè sto sviluppo che mi è uscito mi puzza un pò.

Allora ho preso r(t) e lo sviluppo per entrambe le variabili:
r(t)=k0e−ERT(t)CA(t)

ovvero riscritta per maggiore chiarezza in funzione di x1 e x2

r(t)=k0e−ERx2(t)x1(t)

CODICE:

Codice:

clear all
close all
clc

syms x_1 x_2 a b c d k C_af T_f

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) ==0;

equE= x_1*k*exp(-(b/x_2));

syms R  E

equE= subs(equE, b, E/R)

load('soluzione_x1.mat')
load('soluzione_x2.mat')

T_equ= taylor(equE, [x_1 x_2], [sol_x1 sol_x2])

syms z
equ1= (a)*(C_af-x_1)-z ==0;
equ2= (a)*(T_f-x_2)-(c)*z-(d)*(x_2) ==0;

set_equ=[equ1, equ2];

equZ= subs(set_equ, z, T_equ)

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

Questo è lo sviluppo in serie di taylor credo sia sensato vedo termini esponenziali che rimangono se faccio la derivata parziale in funzione della variabile x1 e in funzione della variabile x2.

CODICE:

Codice:

9.9979155063217651583173638500053*k*exp(-(0.0043336922665508861443061157835321*E)/R) + k*exp(-(0.0043336922665508861443061157835321*E)/R)*(x_1 - 9.9979155063217651583173638500053) - 9.9979155063217651583173638500053*k*exp(-(0.0043336922665508861443061157835321*E)/R)*(x_2 - 230.75011756565803129780229007082)^3*((0.0000000000007642938227379998291357512212758*E^2)/R^2 - (0.00000000035272177890299932004420613223581*E)/R + (0.000018780888661162956802355423359175*E*((0.000000040695295974917565900625149535621*E)/R - (0.000000000058786963150499886674034355372636*E^2)/R^2))/R) - 9.9979155063217651583173638500053*k*exp(-(0.0043336922665508861443061157835321*E)/R)*((0.000000081390591949835131801250299071242*E)/R - (0.00000000017636088945149966002210306611791*E^2)/R^2)*(x_2 - 230.75011756565803129780229007082)^2 - 9.9979155063217651583173638500053*k*exp(-(0.0043336922665508861443061157835321*E)/R)*(x_2 - 230.75011756565803129780229007082)^4*((0.0000000000015285876454759996582715024425516*E)/R - (0.0000000000000033122142289722836809600451822844*E^2)/R^2 + (0.000018780888661162956802355423359175*E*((0.00000000000025476460757933327637858374042527*E^2)/R^2 - (0.00000000017636088945149966002210306611791*E)/R + (0.000018780888661162956802355423359175*E*((0.00000001356509865830585530020838317854*E)/R - (0.000000000014696740787624971668508588843159*E^2)/R^2))/R))/R - (0.000000081390591949835131801250299071242*E*((0.000000040695295974917565900625149535621*E)/R - (0.000000000058786963150499886674034355372636*E^2)/R^2))/R) - 9.9979155063217651583173638500053*k*exp(-(0.0043336922665508861443061157835321*E)/R)*(x_2 - 230.75011756565803129780229007082)^5*((0.000000000000000014354117189256991842248954755283*E^2)/R^2 - (0.0000000000000066244284579445673619200903645688*E)/R - (0.000000081390591949835131801250299071242*E*((0.00000000000025476460757933327637858374042527*E^2)/R^2 - (0.00000000017636088945149966002210306611791*E)/R + (0.000018780888661162956802355423359175*E*((0.00000001356509865830585530020838317854*E)/R - (0.000000000014696740787624971668508588843159*E^2)/R^2))/R))/R + (0.00000000035272177890299932004420613223581*E*((0.000000040695295974917565900625149535621*E)/R - (0.000000000058786963150499886674034355372636*E^2)/R^2))/R + (0.000018780888661162956802355423359175*E*((0.0000000000007642938227379998291357512212758*E)/R - (0.0000000000000011040714096574278936533483940948*E^2)/R^2 + (0.000018780888661162956802355423359175*E*((0.000000000000063691151894833319094645935106317*E^2)/R^2 - (0.000000000058786963150499886674034355372636*E)/R + (0.000018780888661162956802355423359175*E*((0.0000000033912746645764638250520957946351*E)/R - (0.0000000000029393481575249943337017177686318*E^2)/R^2))/R))/R - (0.000000081390591949835131801250299071242*E*((0.00000001356509865830585530020838317854*E)/R - (0.000000000014696740787624971668508588843159*E^2)/R^2))/R))/R) + (0.00018776973796794374142047734691718*E*k*exp(-(0.0043336922665508861443061157835321*E)/R)*(x_2 - 230.75011756565803129780229007082))/R - 1.0*k*exp(-(0.0043336922665508861443061157835321*E)/R)*(x_1 - 9.9979155063217651583173638500053)*(x_2 - 230.75011756565803129780229007082)^4*((0.0000000000015285876454759996582715024425516*E)/R - (0.0000000000000033122142289722836809600451822844*E^2)/R^2 + (0.000018780888661162956802355423359175*E*((0.00000000000025476460757933327637858374042527*E^2)/R^2 - (0.00000000017636088945149966002210306611791*E)/R + (0.000018780888661162956802355423359175*E*((0.00000001356509865830585530020838317854*E)/R - (0.000000000014696740787624971668508588843159*E^2)/R^2))/R))/R - (0.000000081390591949835131801250299071242*E*((0.000000040695295974917565900625149535621*E)/R - (0.000000000058786963150499886674034355372636*E^2)/R^2))/R) - 1.0*k*exp(-(0.0043336922665508861443061157835321*E)/R)*(x_1 - 9.9979155063217651583173638500053)*(x_2 - 230.75011756565803129780229007082)^3*((0.0000000000007642938227379998291357512212758*E^2)/R^2 - (0.00000000035272177890299932004420613223581*E)/R + (0.000018780888661162956802355423359175*E*((0.000000040695295974917565900625149535621*E)/R - (0.000000000058786963150499886674034355372636*E^2)/R^2))/R) - 1.0*k*exp(-(0.0043336922665508861443061157835321*E)/R)*((0.000000081390591949835131801250299071242*E)/R - (0.00000000017636088945149966002210306611791*E^2)/R^2)*(x_1 - 9.997915506321765158317363850[code][/code]0053)*(x_2 - 230.75011756565803129780229007082)^2 + (0.000018780888661162956802355423359175*E*k*exp(-(0.0043336922665508861443061157835321*E)/R)*(x_1 - 9.9979155063217651583173638500053)*(x_2 - 230.75011756565803129780229007082))/R


Secondo te ha senso perchè questo è il risultato che esce soluzione parametrica che è un matrice di 2x5 soluzione rispetto ad x1 è una 1x5 e quella $x_2$ lo stesso. Però è strana perchè la soluzione paramentrica è in funzione di una z che non capisco cos'è. Ne prendo una a campione

Codice:

a*root(- 156250...000*R^5*T_f*a*k - 19899...3*E^5*a*c*k + 1562500.....0000000*R^5*a*k*z + 15625000...0000*R^5*d*k*z - 156250....00000*R^5*T_f*a^2*exp((216....603*E)/(500...0*R)) + 156...00*R^5*a^2*z*exp((216...03*E)/(500...00*R)) - 373....000*E^5*a*c*k*z^2 + 161....000*E^5*a*c*k*z^3 -35094....00*E^5*a*c*k*z^4

ho messo i punti altrimenti era illegibile

[ingres]

Devo guardare a fondo quello che hai postato, ma al momento sono del tuo stesso avviso, ovvero che non sembra granchè convincente.
Dovrebbe trattarsi alla fine dello sviluppo di come variano x1=C_A e x2=T attorno al loro punto di equilibrio in funzione dei parametri, a, b, c, d, k per fissati C_af e Tf , e quindi mi aspetterei, al primo ordine, una formula del tipo:

$Delta x_1= alpha_1*Delta a + beta_1*Delta b + gamma_1*Delta c + delta_1*Delta d+ epsilon1*Delta k$
$Delta x_2= alpha_2*Delta a + beta_2*Delta b + gamma_2*Delta c + delta_2*Delta d+ epsilon2*Delta k$

e non certo una dipendenza polinomiale da z che vedo peraltro definita ad un certo punto come come
(a)*(C_af-x_1)-z ==0
ovvero
z=a*(C_af-x_1)

[MrChopin]

Devo guardare a fondo quello che hai postato, ma al momento sono del tuo stesso avviso, ovvero che non sembra granchè convincente.
Dovrebbe trattarsi alla fine dello sviluppo di come variano x1=C_A e x2=T attorno al loro punto di equilibrio in funzione dei parametri, a, b, c, d, k per fissati C_af e Tf , e quindi mi aspetterei, al primo ordine, una formula del tipo:

$Delta x_1= alpha_1*Delta a + beta_1*Delta b + gamma_1*Delta c + delta_1*Delta d+ epsilon1*Delta k$
$Delta x_2= alpha_2*Delta a + beta_2*Delta b + gamma_2*Delta c + delta_2*Delta d+ epsilon2*Delta k$

e non certo una dipendenza polinomiale da z che vedo peraltro definita ad un certo punto come come
(a)*(C_af-x_1)-z ==0
ovvero
z=a*(C_af-x_1)

no l'ho modificato e l'ho addirittura tolto e mi esce comunque non capisco perchè però non ho fissato Tf e Caf ora ci provo ma non credo cambi molto. Comunque ho fissato tutti i parametri tranne l'esponente b del termine esponenziale, ovvero:

$k_0e^(-E/(Rx_2(t)))x_1(t)=k_0e^(-b/(x_2(t)))x_1(t) $

Codice:

clear all
close all
clc

% Parametri Matlab

F = 1; % Tasso di flusso volumetrico m^3/h
V = 1; % Volume reattore m^3
dH = 5960; % Entalpia 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;
c=dH*k/(rho_cp);
d=UA/(V*rho_cp);

syms x_1 x_2 b

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) ==0;

equE= x_1*k*exp(-(b/x_2));

load('soluzione_x1.mat')
load('soluzione_x2.mat')

T_equ= taylor(equE, [x_1 x_2], [sol_x1 sol_x2])

equ1a= (a)*(C_af-x_1)-T_equ ==0;
equ2a= (a)*(T_f-x_2)-(c)*T_equ-(d)*(x_2) ==0;

set_equ=[equ1a, equ2a];

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

Quello che mi esce è questo

Codice:

(7.68e-413*root(1...0*b^9*z^9 - 8...0*b^8*z^9 + 2...0*b^7*z^9 - 3...0*b^6*z^9 + 2...0*b^5*z^9 - 9...0*b^4*z^9 + 1...0*b^3*z^9 - 1....0*b^2*z^9 - 2...0*b^9*z^8 + 1...0*b^8*z^8 - 5...0*b^7*z^8 + 7...0*b^6*z^8 - 5...0*b^5*z^8 + 2...0*b^4*z^8 - 4...0*b^3*z^8 + 2...0*b^2*z^8 + 2....0*b^9*z^7 - 1...00*b^8*z^7 + 4...0*b^7*z^7 - 7...0*b^6*z^7 + 5...0*b^5*z^7 - 2...0*b^4*z^7 + 4...0*b^3*z^7 - 3...0*b^2*z^7 - 1....0*b^9*z^6 + 9...0*b^8*z^6 - 2...0*b^7*z^6 + 4....0*b^6*z^6 - 3...0*b^5*z^6 + 1....0*b^4*z^6 - 3....0*b^3*z^6 + 2....0*b^5*z^5*exp((2...3*b)/5...0) + 1...0*b^4*z^5*exp((2...3*b)/5....0) - 2....0*b^3*z^5*exp((2...3*b)/5...0) + 1....0*b^2*z^5*exp((2...3*b)/5....0) - 1....0*b*z^5*exp((2...3*b)/5...0) + 3...0*b^8*z^5 + 1....0*b^7*z^5 - 1....0*b^6*z^5 + 1....0*b^5*z^5 - 6...0*b^4*z^5 + 1...0*b^3*z^5 - 1...0*b^2*z^5 + 7...0*b*z^5 + 4....0*b^5*z^4*exp((2...3*b)/5...0) - 2...0*b^4*z^4*exp((2...0) + 3...0*b^3*z^4*exp((2...0) - 1...0*b^2*z^4*exp((2...3*b)/5...0) + 2...0*b*z^4*exp((2...3*b)/5...0) - 9...0*b^9*z^4 + 7...0*b^8*z^4 - 2...0*b^7*z^4 + 4....0*b^6*z^4 - 3...0*b^5*z^4 + 1....0*b^4*z^4 - 3....0*b^3*z^4 + 3...0*b^2*z^4 -... Output truncated. Text exceeds maximum line length for Command Window display.

E penso che continui all'infinito non c'è un modo per renderla più compatta e leggibile? Tipo una riduzione ai minimi termini mettendo in evidenza tutto ciò che è possibile?

[ingres]

Devo vedere con Matlab per bene lo script che mi hai allegato, anche se una soluzione che parte con 7.68e-413 non mi sembra che possa funzionare.

Comunque, per altra via, ho calcolato ad es. i coefficienti del primo ordine di influenza sulla soluzione di equilibrio per a=F/V che dovrebbero essere:

$alpha_1 = 0.33$
$alpha_2 = 3.42$

Guarda se riesci a riottenere questo risultato.

[MrChopin]

Devo vedere con Matlab per bene lo script che mi hai allegato, anche se una soluzione che parte con 7.68e-413 non mi sembra che possa funzionare.

Comunque, per altra via, ho calcolato ad es. i coefficienti del primo ordine di influenza sulla soluzione di equilibrio per a=F/V che dovrebbero essere:

$alpha_1 = 0.33$
$alpha_2 = 3.42$

Guarda se riesci a riottenere questo risultato.

Aspetta non ho capito come li hai ottenuti

[ingres]

Semplicissimo. Ho variato F/V di una quantità $Delta$ e ho trovato i nuovi valori di equilibrio. Quindi ho fatto il rapporto tra la variazione del punto di equilibrio e il $Delta$.

[MrChopin]

Semplicissimo. Ho variato F/V di una quantità $Delta$ e ho trovato i nuovi valori di equilibrio. Quindi ho fatto il rapporto tra la variazione del punto di equilibrio e il $Delta$.

Ah ok alla fine stai facendo un sorta algoritmo evoluzionistico patter serch. Effittivamente sarebbe interessante creare uno codice patter search per creare un dominio di variazioni parametriche che non mi fanno scomparire il punto di equilibrio. Creando una sorta di dominio sicuro una sorta di LAS euristica (Local asymptotic stability). Sempre se non ho detto fesserie e non ha alcun senso quello che ho pensato. Però io devo cercare la formula parametrica "bruta" della serie $x=a/b$.