Simulacion Sistemas Biologicos en Matlab

UNIVERSIDAD AUTONOMA CHAPINGO SIMULACION DE SISTEMAS BIOLOGICOS PRIMERA TAREA

1. Programe el siguiente modelo matemático en Matlab-Simulink dx1 dt dx2 dt

 ax1  bu1

 cx1  dx2  eu2

Donde a  2 , b  5 , c  3 , d  6 , e  2 . valor de 0 a 2 en el tiempo

t  3 s.

u1 (t )

es una función escalón que cambia de un

u2 (t ) es una función senoidal con una frecuencia de 0.1

hz. Las condiciones iniciales son x1 (0)  0 , x2 (0)  0 . Utilice el método de integración de Runge Kutta cuarto orden (ode45).

Primeramente se elabora el modelo en simulink (ver “problema1sim.mdl ”) o Figura 1.

Figura 1. Modelo en simulink A continuación se presenta el programa hecho en matlab o ver “problema1.m ”: clear all a=-2; b=5; c=3; d=-6; e=2; x1=0; x2=0; opciones=simset('solver' opciones=simset( 'solver', ,'ode45' 'ode45', ,'reltol' 'reltol',1e-8); ,1e-8); [t,x,y]=sim('problema1sim' [t,x,y]=sim( 'problema1sim',[0 ,[0 300],opciones); subplot(2,1,1); plot(t,x(:,1)); grid

PRESENTA. ING. PEDRO BARRERA PUGA

UNIVERSIDAD AUTONOMA CHAPINGO SIMULACION DE SISTEMAS BIOLOGICOS PRIMERA TAREA xlabel('Tiempo xlabel('Tiempo en segundos') segundos') ylabel('Y' ylabel('Y') ) subplot(2,1,2); plot(t,x(:,2)); grid xlabel('Tiempo xlabel('Tiempo en segundos') segundos') ylabel('Y' ylabel('Y') )

Una vez terminado se corre la simulación obteniéndose los resultados de la Figura 1.2 :

Figura 1.2 Simulación Simulación de las dos ecuaciones ecuaciones diferenciales. 2. Programe en Matlab-Simulink el siguiente modelo simplificado para el crecimiento de biomasa dw  c1 f (T , I )(1  e c2w )  c3w 2 dt donde w es peso de la planta en kg, c1  0.33 , c2  0.01 , c3  0.4 , w(0) (0)  0.00 0.0002 02 , todas son constantes en una escala de hora, f ( T, I ) es una función que hace depender la tasa de crecimiento de la biomasa en función de la temperatura y la radiación. Simule durante una semana y durante 75 días. Considere que f ( T, I )  1 .

Modelo en simulink simulink se muestra en la Figura 2.1 o ver “problema2.mdl ”

Figura 2.1 Modelo en simulink del problema 2 PRESENTA. ING. PEDRO BARRERA PUGA

UNIVERSIDAD AUTONOMA CHAPINGO SIMULACION DE SISTEMAS BIOLOGICOS PRIMERA TAREA Programa en Matlab, ver “problema2a.m ”. clear % PROBLEMA 2. c1=0.33; c2=0.01; c3=0.4; d=24*7; % simulacion de una semana w0=0.0002; opciones=simset('solver' opciones=simset( 'solver', ,'ode45' 'ode45', ,'reltol' 'reltol',1e-8); ,1e-8); [t,x,y]=sim('problema2' [t,x,y]=sim( 'problema2',[0 ,[0 d],opciones); subplot(2,1,1); plot(t/24,x(:,1),'r' plot(t/24,x(:,1), 'r'); ); xlabel('Tiempo xlabel('Tiempo en días'); días'); ylabel('Masa ylabel('Masa en kg'); kg'); title('CRECIMIENTO title('CRECIMIENTO DE BIOMASA DE UNA SEMANA'); SEMANA' ); grid d1=24*75; %simulacion de 75 dias. opciones=simset('solver' opciones=simset( 'solver', ,'ode45' 'ode45', ,'reltol' 'reltol',1e-8); ,1e-8); [t,x,y]=sim('problema2' [t,x,y]=sim( 'problema2',[0 ,[0 d1],opciones); subplot(2,1,2); plot(t/24,y,'r' plot(t/24,y, 'r'); ); xlabel('Tiempo xlabel('Tiempo en días'); días'); ylabel('Masa ylabel('Masa en kg'); kg'); title('CRECIMIENTO title('CRECIMIENTO DE BIOMASA DE 75 DIAS'); DIAS' ); grid

Resultados de la simulación se muestran en la Figura 2.2:

Figura 2.2 Comportamiento de la biomasa durante una semana y 75 días.


UNIVERSIDAD AUTONOMA CHAPINGO SIMULACION DE SISTEMAS BIOLOGICOS PRIMERA TAREA 3. Programe el siguiente modelo matemático en Matlab-Simulink dT a dt



1 C gr

 R  kco (Ta  To )  kvent (Ta  To )  H  -1

C gr : Capacidad de calor total del invernadero (32000 J °C ) k co

: Coeficiente de transferencia de calor de la cubierta (7.0 J s

-1

-1

°C )

k vent : Coeficiente de transferencia de calor debido a la ventilación (15.0 J s

H R T o

-1

°C-1)

-1

: Calor producido por el sistema de calefacción (J s ) : Radiación solar (200 J s -1) : Temperatura del aire fuera del invernadero (5°C)

T a (0)  0 °C

a) Transforme la escala del modelo de segundos a horas. b) Determine la solución en estado estado estacionario analíticamente si no hay calefacción (H=0) c) Verifique la solución analítica numéricamente usando Matlab d) Utilice la solución en estado estacionario como condición inicial y suponga que -1 el sistema de calefacción se activa (H=200 J s ) e) Determine el nuevo estado estacionario. f) Realice simulaciones y genere graficas de los resultados.

El modelo en simulink se observa en la Figura 3.1 º ver “problem3.mdl ”:

Figura 3.1 Simulación de la Temperatura del aire Programa en Matlab o ver “problema3.m ”: %PROBLEMA 3 (a) clear Cgr=32000; Kc=7.0; Kv=15.0; H=200; R=100; To=5; Ta=0; pedro=simset('solver' pedro=simset( 'solver', ,'ode23' 'ode23', ,'reltol' 'reltol',1e-6); ,1e-6); [t,x,y]=sim('problem3' [t,x,y]=sim( 'problem3',[0 ,[0 60*60*8],pedro); subplot(2,1,1); plot(t/(60*60),x(:,1),'r' plot(t/(60*60),x(:,1), 'r'); ); grid xlabel('tiempo xlabel('tiempo en horas') horas') ylabel('temperatura ylabel('temperatura en ºC') ºC')


UNIVERSIDAD AUTONOMA CHAPINGO SIMULACION DE SISTEMAS BIOLOGICOS PRIMERA TAREA title('TEMPERATURA title('TEMPERATURA DENTRO DEL INVERNADERO H=200 Js-1 ') ') hold on on; ; H=0; pedro=simset('solver' pedro=simset( 'solver', ,'ode23' 'ode23', ,'reltol' 'reltol',1e-6); ,1e-6); [t,x,y]=sim('problem3' [t,x,y]=sim( 'problem3',[0 ,[0 60*60*8],pedro); subplot(2,1,2); plot(t/(60*60),y,'r' plot(t/(60*60),y, 'r'); ); grid xlabel('tiempo xlabel('tiempo en horas') horas') ylabel('temperatura ylabel('temperatura en º C') C') title('TEMPERATURA title('TEMPERATURA DENTRO DEL INVERNADERO H=0 Js-1' Js-1') )

los resultados de la Simulación para una H=200 y H=0 se muestran en la Figura 3.2:

Figura 3.2 Comportamiento de la temperatura del aire con dos valores de Humedad En estado estacionario dT a  0; dt

1 C gr

 R  k CO T a  T o   k vent T a  T o   0

T a  T o 

 R , Ta es igual a:  k co  k ven


UNIVERSIDAD AUTONOMA CHAPINGO SIMULACION DE SISTEMAS BIOLOGICOS PRIMERA TAREA T a 

 R  T o , Para tener el valor de Ta se sustituyen los valores de las  k co  k ven

constantes.

 200 * 3600 T a   5  14.09  (7 * 3600)  (15 * 3600)

J  s 1 J  s 1 º C 1

º C

b) Verifique la solución solución analítica numéricamente usando usando Matlab % Se empleo “trim” para realizar este procedimiento H=0 estacionario=trim( 'problem3' 'problem3')) % c)se uso "trim". H=0 eet=trim('problem3' eet=trim('problem3') ) % d) Ta(0)= eet H=200; Ta=eet; eet1=trim('problem3' eet1=trim('problem3') )

% e) H=200 H=200 Ta=23.18; pedro=simset('solver' pedro=simset( 'solver', ,'ode23' 'ode23', ,'reltol' 'reltol',1e-6); ,1e-6); [t,x,y]=sim('problem3' [t,x,y]=sim( 'problem3',[0 ,[0 3600*dia],pedro); plot(t/(3600),y(:,1)) grid xlabel('tiempo xlabel('tiempo [h]') [h]') xlabel('tiempo xlabel('tiempo [horas]') [horas]') ylabel('Temp ylabel('Temp de nuevo E. Estacionario [ºC]') [ºC]' )

La Figura 3.3 presenta el comportamiento de la temperatura con el nuevo estado estacionario calculado.

Figura 3.3 Nuevo estado estacionario de la temperatura. PRESENTA. ING. PEDRO BARRERA PUGA

UNIVERSIDAD AUTONOMA CHAPINGO SIMULACION DE SISTEMAS BIOLOGICOS PRIMERA TAREA 4. Programe el siguiente modelo matemático en Matlab-Simulink dT a dt dT s



dt C s



1 C gr

1 C s

 R  kco (Ta  To )  kvent (Ta  To )  kas (Ta  Ts )  H 

kas (Ta  Ts )  ksd (Ts  Tsd ) 

: Capacidad de calor de la capa de suelo (120000 J)

k as :

Coeficiente de transferencia de calor entre el suelo y el aire (5.0 J s -1°C-1 )

k sd :

Coeficiente de transferencia de calor entre la capa superior y capa profunda del

suelo (2.0 J s-1°C-1 ). T sd : Temperatura de la capa profunda del suelo (10 °C) a) Cambie la escala del modelo de segundos a horas b) Determine la solución en estado estacionario del modelo para H=0. % (b) H=0 Ta_1=trim('problesim4c' Ta_1=trim('problesim4c') )

c) Aplique un valor de H=200 mediante una señal escalón en el tiempo t  1 hora. Simule durante varias horas el comportamiento de las dos variables. El modelo en simulink se presenta en la figura 4.1 o ver también “problesim4c.mdl ”:

Figura 4.1 Modelo de temperatura del aire y temperatura del suelo Programa en matlab o “problema4c.m ”: %PROBLEMA 4 (c) clear Cgr=32000; Kco=7.0; Kve=15.0; H=200; %Valor que se cambia en la función escalón H=0.


UNIVERSIDAD AUTONOMA CHAPINGO SIMULACION DE SISTEMAS BIOLOGICOS PRIMERA TAREA Rad=200; To=5; Cs=120000; Kas=5.0; Ksd=2.0; Tsd=10.0; Ta=0; Ts=0; opciones=simset('solver' opciones=simset( 'solver', ,'ode45' 'ode45', ,'reltol' 'reltol',1e-8); ,1e-8); [t,x,y]=sim('problesim4c' [t,x,y]=sim( 'problesim4c',[0 ,[0 60*60*50],opciones); %Aquí se cambia la escala de segundos a horas subplot(2,1,1); plot(t/(60*60),y(:,1)) grid title('TEMPERATURA title('TEMPERATURA DEL SUELO') SUELO' ) xlabel('Tiempo xlabel('Tiempo en horas') horas') ylabel('Temp. ylabel('Temp. [Celsius]') [Celsius]') subplot(2,1,2) plot(t/(60*60),y(:,2),'r' plot(t/(60*60),y(:,2), 'r') ) grid title('TEMPERATURA title('TEMPERATURA DEL AIRE') AIRE') xlabel('Tiempo xlabel('Tiempo en horas') horas') ylabel('Temp. ylabel('Temp. [Celsius]') [Celsius]')

En la figura 4.2 se presenta los graficos de la simulacion de temperatura del aire y temperatura del suelo:

Figura 4.2 Comportamiento de la temperatura del aire y temperatura del suelo en 50 días


UNIVERSIDAD AUTONOMA CHAPINGO SIMULACION DE SISTEMAS BIOLOGICOS PRIMERA TAREA 5. Genere una función para la radiación solar que cambie de acuerdo con el tiempo del día, como sigue: si t  18:00 y t  6:00 horas R  0 R R0 sin(2 t) si t  6 y t  18 Con R0  200 (J s-1) Verifique el comportamiento de su programa. Combine la función para la radiación solar con el modelo generado en el problema -1 4 y 5. Con H=200 J s , realice una simulación durante varios días. Genere graficas de los resultados. Primero se generó la función para la temperatura, t emperatura, la cual fue utilizada en la función “ Math function” o ver “funt5.m”: function R=funct5(u) % nombre da la función para utilizar con MATH FUNCTION de Simulink t=u(1); Rad=u(2); if ((t>0.75) | (t<0.25)) R=0; else R=Rad*sin(pi*t); end

Se probó la función para un día de acuerdo al siguiente programa de matlab o ver “ problema5a.m ”: clear Rad=200; %rdh=3600; dia=1 t=[0:1/24:1]; opciones=simset('solver' opciones=simset( 'solver', ,'ode45' 'ode45', ,'reltol' 'reltol',1e-6); ,1e-6); [t,x,y]=sim('prue5a' [t,x,y]=sim( 'prue5a', , t,opciones); subplot(2,1,1); plot(t,y); grid; title('VARIACION title('VARIACION DE LA RADIACIÓN') RADIACIÓN') xlabel('Tiempo xlabel('Tiempo en días') días') ylabel('Radicación ylabel('Radicación Js-1') Js-1')

El modelo de simulink (Figura 5.1) que se llama a matlab es “ prue5a.mdl ”.

Figura 5.1 Diagrama para la función de radiación solar PRESENTA. ING. PEDRO BARRERA PUGA

UNIVERSIDAD AUTONOMA CHAPINGO SIMULACION DE SISTEMAS BIOLOGICOS PRIMERA TAREA Luego en la Figura 5.2 se genero la combinación de la función de radiación solar con el modelo generado en el problema 4 y se hace una simulación durante varios días con H=200 J s -1 El modelo en simulink es el siguiente o ver “ problesim5b.mdl ”:

Figura 5.2 Diagrama donde se combina combina radiación solar , temperatura del aire y del suelo El programa generado en matlab (ver “problema5b.m” en combinación se muestra a continuación: %PROBLEMA 5 (b) clear hr=3600; Cgr=32000; Kco=7.0*hr; Kve=15.0*hr; H=200*hr; Rad=200*hr; To=5; Cs=120000; Kas=5.0*hr; Ksd=2.0*hr; Tsd=10.0; Ta=0; Ts=0; %dia=1; opciones=simset('solver' opciones=simset( 'solver', ,'ode45' 'ode45', ,'reltol' 'reltol',1e-6); ,1e-6); [t,x,y]=sim('problesim5b' [t,x,y]=sim( 'problesim5b',[0 ,[0 15],opciones); subplot(2,1,1); plot(t,y(:,1)) title('Radiacion title('Radiacion Solar') Solar') xlabel('tiempo xlabel('tiempo [DIAS]') [DIAS]') ylabel('Radiacion ylabel('Radiacion [Js-1]') [Js-1]') grid subplot(2,1,2);


UNIVERSIDAD AUTONOMA CHAPINGO SIMULACION DE SISTEMAS BIOLOGICOS PRIMERA TAREA plot(t,y(:,2)) title('Variacion title('Variacion de Temperatura') Temperatura') xlabel('tiempo xlabel('tiempo [DIAS]') [DIAS]') ylabel('temperatura ylabel('temperatura [Grados Centígrados]') Centígrados]') grid

Y los resultados de la simulación para 15 días se presentan en la Figura 5.3:

Figura 5.3 Comportamiento de la temperatura y radiación por 15 días 6. Modifique el problema 2 considerando ahora la función f ( T, I ) como sigue I , con f (T , I )  c4 f (Ta ) c5  I f (T )  0 si T a  10 f( T) 

1 15

( Ta  10) si 10  T a  25

1

( Ta  25) si 25  T a  30 5 f (T )  0 si T a  30 f( T)  1 

Donde la temperatura está dada en °C, c4  4.5 , c5  250000 Simule el sistema durante una semana y también durante 75 días usando valores variables de la temperatura.

Se elaboró el programa en Matlab-Simulink quedando como se ilustra en la Figura 6.1 o ver pro66.mdl ”:



Figura 6.1 Modelo de biomasa biomasa y temperatura Como se observa, en la Figura 6.1. Se utilizaron dos funciones que fueron elaboradas en otro archivo, la función de Temperatura y la de Radiación solar, las cuales son: biom1.m para biomasa, y funT.m para temperatura.

Y la función de Temperatura:


UNIVERSIDAD AUTONOMA CHAPINGO SIMULACION DE SISTEMAS BIOLOGICOS PRIMERA TAREA A continuación estas dos funciones son llamadas por el bloque de simulink “ Matlab Math Function”.

Se elaboró el script en Matlab quedando de la siguiente manera o ver “problema6.m ”: %SOLUCION AL PROBLEMA 6% clear c1=0.33; c2=0.01; c3=0.4; c4=4.5; fa=3600; c5=250000; I=200*fa; w0=0.0002; day1=24*7; Ta=20; Dia=1; Tmaxima=30; Tminima=15; opciones=simset('solver' opciones=simset( 'solver', ,'ode45' 'ode45', ,'reltol' 'reltol',1e-8); ,1e-8); [t,x,y]=sim('pro66' [t,x,y]=sim( 'pro66',[0 ,[0 day1],opciones); subplot(2,1,1) plot(t/24,y(:,1),'b' plot(t/24,y(:,1), 'b') ) title('Biomasa title('Biomasa de acuerdo a la temperatura [7 días]') días]' ) xlabel('tiempo xlabel('tiempo [días]') [días]') ylabel('biomasa' ylabel('biomasa') ) grid day=24*75 opciones=simset('solver' opciones=simset( 'solver', ,'ode45' 'ode45', ,'reltol' 'reltol',1e-6); ,1e-6); [t,x,y]=sim('pro66' [t,x,y]=sim( 'pro66',[0 ,[0 day],opciones); subplot(2,1,2) plot(t/24,y(:,1)) title('Biomasa title('Biomasa de acuerdo a la temperatura [75 días]') días]' ) xlabel('tiempo[h]' xlabel('tiempo[h]') ) ylabel('biomasa ylabel('biomasa ') ') grid

La Figura 6.2 Muestra la simulación:



Figura 6.2 Comportamiento de la biomasa por 75 días y una semana 7. Combine el modelo de crecimiento de biomasa (obtenido en el problema 6) con el modelo de la temperatura del aire del problema 4. Este es un ejemplo de un biosistema rígido que requiere de métodos especiales de integración. Utilice los métodos de integración para sistemas rígidos (stiff) disponibles en Matlab. Simule el comportamiento del sistema completo durante 75 días. Utilice un valor de H  350 J/s.

De la misma manera, se elaboro la combinación del problema 6 y problema 5, y se utilizo la misma función “biom1.m”, quedando el programa en simulink y ver “problesim7.mdl” como se muestra en la figura 7.1 .

Figura 7.1, Programa de bloques de combinación de biomasa y temperatura. PRESENTA. ING. PEDRO BARRERA PUGA

UNIVERSIDAD AUTONOMA CHAPINGO SIMULACION DE SISTEMAS BIOLOGICOS PRIMERA TAREA A continuación se elaboro el script en matlab, ver también “problema7.m ”: %SOLUCION AL PROBLEMA 7% clear fa=3600; Cgr=32000; Kco=7.0*fa; Kve=15.0*fa; H=350*fa; Rad=200*fa; To=5; Cs=120000; Kas=5.0*fa; Ksd=2.0*fa; Tsd=10.0; Ta=0; Ts=0; c1=0.33; c2=0.01; c3=0.4; c4=4.5; c5=250000; I=200*fa; w0=0.0002; day1=24*200; Ta=20; Dia=1; Tmaxima=30; Tminima=15; opciones=simset('solver' opciones=simset( 'solver', ,'ode23' 'ode23', ,'reltol' 'reltol',1e-8); ,1e-8); [t,x,y]=sim('problesim7' [t,x,y]=sim( 'problesim7',[0 ,[0 day1],opciones); plot(t/(24),y(:,3),'b' plot(t/(24),y(:,3), 'b') ) title('Biomasa title('Biomasa de acuerdo a la temperatura [200 días]') días]' ) xlabel('tiempo xlabel('tiempo [días]') [días]') ylabel('biomasa' ylabel('biomasa') ) grid

Y Se obtuvieron los resultados de la simulación en la Figura 7.2:



Figura 7.2 Comportamiento de la biomasa, durante 200 días La simulación se elaboró para 200 días, ya que para 75 días no se observa claramente el comportamiento de la biomasa. 8. Combine el modelo de crecimiento de biomasa (obtenido en el problema 6) con el modelo de la temperatura del aire del problema 5. Este es un ejemplo de un biosistema rígido que requiere de métodos especiales de integración. Utilice los métodos de integración para sistemas rígidos (stiff) disponibles en Matlab. Simule el comportamiento del sistema completo durante 75 días. Utilice H  350 J/s y otros valores.

Posteriormente se programo el programa en simulink el cual se muestra en la Figura 8.1 o ver también “problesim8.mdl”



Figura 8.1 Modelo de combinación de modelo de biomasa y de temperatura del aire. De igual manera se utilizo la función de radiación y de temperatura “funtrad.m” la cual es:

Después se elaboro el script en Matlab( ver “problema8.m ”: %SOLUCION AL PROBLEMA 8% clear hr=3600; Cgr=32000; Kco=7.0*hr; Kve=15.0*hr; H=350*hr; Rad=200*hr; To=5; Cs=120000; Kas=5.0*hr; Ksd=2.0*hr; Tsd=10.0; Ta=0; Ts=0;


UNIVERSIDAD AUTONOMA CHAPINGO SIMULACION DE SISTEMAS BIOLOGICOS PRIMERA TAREA c1=0.33; c2=0.01; c3=0.4; c4=4.5; fa=3600; c5=250000; I=200*fa; w0=0.0002; day1=24*75; Ta=20; opciones=simset('solver' opciones=simset( 'solver', ,'ode45' 'ode45', ,'reltol' 'reltol',1e-8); ,1e-8); [t,x,y]=sim('problesim8' [t,x,y]=sim( 'problesim8',[0 ,[0 day1],opciones); plot(t/24,y(:,3)) title('Biomasa title('Biomasa de acuerdo a la temperatura [75 días]') días]' ) xlabel('tiempo[días]' xlabel('tiempo[días]') ) ylabel('biomasa ylabel('biomasa ') ') grid

A continuación en la Figura 8.2 Se muestra los resultados de la la simulación de la combinación del modelo.

Figura 8.2 Simulación del comportamiento de crecimiento de biomasa y temperatura 9. Linealización a. Linealice el modelo del problema 5 en forma analítica alrededor de la solución estado estacionario para H=0. b. Obtenga tanto la solución en estado estacionario cuando, H=0, usando MATLAB. c. Compare las soluciones.



El programa en Matlab se muestra a continuación o ver también “problema9.m ”: clear hor=3600; Cgr=32000; Kco=7.0*hor; Kve=15.0*hor; H=350*hor; Rad=200*hor; To=5; Cs=120000; Kas=5.0*hor; Ksd=2.0*hor;


UNIVERSIDAD AUTONOMA CHAPINGO SIMULACION DE SISTEMAS BIOLOGICOS PRIMERA TAREA Tsd=10.0; Ta=0; Ts=0; H=350*hor [A,B,C,D]=linmod('problesim4c' [A,B,C,D]=linmod( 'problesim4c') ) B=[1;0]; D=[0;0];

En este problema únicamente lo que se hizo fue cambiar la H, por H=350, y agregar los componentes de la matriz. matr iz. Obteniéndose los siguientes resultados:

10. Programe en Matlab-Simulink el modelo simplificado para crecimiento de tomate descrito en la sección VIII del capítulo Simulation of biological processes por James W. Jones and Joep C. Luyten

Primeramente se elaboro el programa en simulink para determinar el número de nodos por planta, y como como el número de nodos por planta depende depende de la temperatura entonces se programo una función función para las condiciones de temperatura temperatura que se nos indico y para interpolar los valores intermedios se utilizo la función interp1 (Ver función “temperatura.m”): dN dt

 r m r (T ) (26)

Donde: dN dt

:

Número de nodos con respecto al tiempo.

Tasa de aparición de hojas ho jas r (T ) : Es una función de la temperatura. r m

:

Las condiciones para elaborar la función son:


UNIVERSIDAD AUTONOMA CHAPINGO SIMULACION DE SISTEMAS BIOLOGICOS PRIMERA TAREA Es cero cuando la Temperatura sea menor de 8°C o mayor de 50°C y tiene los valores de 0.55, 1.0 y 1.0 cuando la temperatura es 12, 30 y 35 °C respectivamente. Para el ejemplo se dice que la tasa de aparición aparición de hojas es máxima entre 30 y 35 °C. r (T ) :

Posteriormente se elaboro la la determinación del índice índice de área foliar con la la ecuación 32: L   ( /  ) ln1  exp exp  ( N  nb )(32)

Donde: L: Índice de área foliar (m 2 de hoja) /(m 2 de suelo) 2  : Densidad de planta, numero/m N: Número de hoja.  ,  ,  y nb : Son coeficientes empíricos de la ecuación expolineal   delta  0.074   beta  0.38 nb  nb  13.3

Se elaboro otra función ( ver “densidadluz.m ”) para determinar la densidad de flujo

luminoso en lo alto del dosel, empleando para ello la ecuación 31. I 0  I m sin2 (t h  6) / 24 …….(31)


UNIVERSIDAD AUTONOMA CHAPINGO SIMULACION DE SISTEMAS BIOLOGICOS PRIMERA TAREA Donde: t h :

Es el tiempo solar en horas. I : Máxima densidad de flujo luminoso (en la noche) 2 I 0 : Densidad del flujo luminoso en lo alto del dosel, (µmol de fotones)/[(m de suelo)*s] m

Asumiendo 12 horas diarias se tiene que: 6  t h  18 A continuación continuación se presenta la función para la densidad de flujo flujo luminoso ( programa “densidadluz.m ”):

Después se calculo el factor de reducción de fotosíntesis con la ecuación 30 p(T )  1  (( h  T ) /( h   1 ))2 ……….(30)

Donde: p(T ) : Es el factor de reducción de fotosíntesis adimencional.

 h  fih  30.0 C : Es la temperatura a la cual la fotosíntesis de la hoja es máxima.  1  fi1  5.0 C : Es la temperatura bajo la cual la fotosíntesis de la hoja es cero. T: temperatura en °C que varía entre 6 y 18 °C. ( Ver función “temperatura2.m ”)


UNIVERSIDAD AUTONOMA CHAPINGO SIMULACION DE SISTEMAS BIOLOGICOS PRIMERA TAREA Con las ecuaciones 30, 31 y 32, se prosiguió prosiguió a elaborar la Tasa de fotosíntesis de la ecuación 29. Pg  D

 Cp(T ) K



  ………(29)  exp( ) ( 1 )  KI  KL   m C 0  

ln 

 KI 0  (1  m) C

Donde: D=0.108= coeficiente para convertir la fotosíntesis de (µmol de CO2)/(m2*s) a (g CH2O)/(m2*h) 2  =tau=0.0664: conductancia del CO2 en la hoja, (µmol de CO 2)/[(m de hoja)*s] C=350: Concentración de CO2 del aire, (µmol (µ mol de CO2)/(mol de aire) (µ mol de  =alfa=0.056: Eficiencia de utilización de la luz por la hoja, (µmol de CO 2)/ (µmol fotones) K=0.58: Coeficiente de extinción de luz en el dosel, adimensional. m  0.10 : Coeficiente de transmisión de luz de las hojas, adimensional. L: Índice de área foliar del dosel, (m2 de hoja ) /(m2 de suelo) Se programo la tasa de respiración de mantenimiento R m con la ecuación 27. Rm  k m exp exp 0.0693(T  25)………(27)

Donde: Rm: Tasa de respiración de mantenimiento, (g CH 2O)/[(g de tejido)*h] T: Temperatura en °C k m  km  0.0006 : Tasa de respiración a 25 °C, (g CH 2O)/[(g de tejido)*h] Se elaboro la programación en simulink de la tasa de crecimiento de materia seca de la ecuación 28. dW dt

 E ( Pg  RmW ) …………..(28)

Donde: dW dt

: Es la tasa de crecimiento de la cosecha del peso seco, ( g de tejido)/(m 2*h)

W: total de materia seca por planta, g/m2 E=0.70: eficiencia de conversión de CH2O a tejido de planta, ( g de tejido)/(g de CH 2O) Pg: Tasa de fotosíntesis bruta del dosel, (g de CH 2O )/(m2*h)


UNIVERSIDAD AUTONOMA CHAPINGO SIMULACION DE SISTEMAS BIOLOGICOS PRIMERA TAREA Finalmente se programo la la ecuación 33, para simular simular la biomasa del dosel y la biomasa biomasa de la raíz por separado. dW c dt

 E ( Pg  RmW ) f c ( N ),

dW r dt

 E ( Pg  RmW )[1  f c ( N )] ……………. (33)

Donde: W: total de materia seca por planta, g/m2 E=0.70: eficiencia de conversión de CH2O a tejido de planta, ( g de tejido)/(g de CH 2O) Pg: Tasa de fotosíntesis bruta del dosel, (g de CH 2O )/(m2*h) Rm: Tasa de respiración de mantenimiento, (g CH 2O)/[(g de tejido)*h] f c ( N ) =0.85 dW c dt dW r dt

: Es la tasa de crecimiento de la cosecha del peso seco del dosel, ( g de tejido)/(m 2*h) : Es la tasa de crecimiento de la cosecha del peso seco de la raíz, ( g de tejido)/(m 2*h)

Enseguida se muestra el script elaborado en Matlab y ver “preoblema10.m ”. %SOLUCION DEL PROBLEMA 10 clear dN=0.0; E=0.70; fc=0.85; rm=0.021*24; Dia=24; dWc=0; dWr=0; dW=0; Tmaxima=30; Tminima=20; fih=30.0; fi1=5.0; D=0.108*24; tau=0.0664; delta=0.074; beta=0.38; nb=13.3; km=0.0006*24; C=350; K=0.58; alfa=0.056; m=0.10; Im=1200; ro=4.0;


UNIVERSIDAD AUTONOMA CHAPINGO SIMULACION DE SISTEMAS BIOLOGICOS PRIMERA TAREA t=[0:1/24:80]; % AQUI SE SIMULO PARA 80 DÍAS COMO INDICA EL ARTÍCULO pedro=simset('solver' pedro=simset( 'solver', ,'ode23' 'ode23', ,'reltol' 'reltol',1e-8); ,1e-8); [t,x,y]=sim('problem10sim' [t,x,y]=sim( 'problem10sim',t ,t ,pedro); figure(1); plot(t,y(:,2)); xlabel('tiempo xlabel('tiempo [dias]'); [dias]'); ylabel('Indice ylabel('Indice de Area Foliar[(m2 de hoja)(m2 de suelo)]'); suelo)]' ); grid title('IAF' title('IAF'); ); figure(2); plot(t,y(:,5)); grid xlabel('tiempo xlabel('tiempo [dias]'); [dias]'); ylabel('Tasa ylabel('Tasa de fotosintesis [(g de CH2O )/(m2)]'); )/(m2)]' ); figure(3); plot(t,y(:,7)); grid xlabel('tiempo xlabel('tiempo [dias]'); [dias]'); ylabel('biomasa ylabel('biomasa de la raiz [g de hoja/ m2 de hoja]'); hoja]' ); figure(4); plot(t,y(:,8)); grid xlabel('tiempo xlabel('tiempo [dias]'); [dias]'); ylabel('biomasa ylabel('biomasa del dosel [g de hoja/ m2 de hoja]'); hoja]' ); figure(5); plot(t,y(:,9)); grid xlabel('tiempo xlabel('tiempo [dias]'); [dias]'); ylabel('biomasa ylabel('biomasa total [g de hoja/ m2 de hoja]'); hoja]' ); Finalmente presenta el modelo en simulink ( ver “problema10sim”).



Se obtuvieron las siguientes figuras de la simulación:






Simulacion Sistemas Biologicos en Matlab

Recommend Documents