Problemas de control en procesos de eutrofización

XX Congreso de Ecuaciones Diferenciales y Aplicaciones X Congreso de Matem´ atica Aplicada Sevilla, 24-28 septiembre 2007 (pp. 1–8) Problemas de control en procesos de eutrofizaci´on L. J. ´ Alvarez-V´ azquez 1, F. Javier Fern´ andez2, R. Mu˜ noz-Sola3 1Dpto. Mat. Apl. II. E.T.S.I. Telecomunicaci´on. Universidad de Vigo. E-mail: [email protected]. 2Dpto. Mat. Apl. Universidad de Santiago de Compostela. E-mail: [email protected]. 3Dpto. Mat. Apl. Universidad de Santiago de Compostela. E-mail: [email protected]. Palabras clave: Control ´ Optimo, Eutrofizaci´on, Biorreactores Resumen El aumento de los niveles de nitr´ogeno y f´osforo (nutrientes principales para determinadas algas) en grandes masas de agua debido a la acci´on contaminante del hombre, puede traer consigo una disminuci´on de la calidad del agua (procesos de eutrofizaci´on). En este trabajo estudiamos un m´etodo de control de esta clase de procesos, basado en el uso de biorreactores (tanques en donde se desarrollar´an una serie de procesos que reducir´an los niveles de nitr´ogeno y/o f´osforo del agua por debajo de unos umbrales m´aximos). Matem´aticamente estos m´etodos de control se pueden formular como problemas de control ´optimo con restricciones sobre el estado en donde las ecuaciones de estado son un tipo de sistema de ecuaciones cuasilineales. 1. Las ecuaciones de estado. Un sistema de ecuaciones cuasilineales. Las ecuaciones que modelan los procesos de eutrofizaci´on se pueden considerar como un caso particular del siguiente sistema de ecuaciones no lineal:                  ∂ul ∂t −Pn i=1 ∂ ∂xial i(t, x, ul,∇ul) + cl(t, x, u,∇u) = glen Q, ν·al(t, x, ul,∇ul) + bl(t, x, ul) = hlsobre Σ, ul(0) = ul 0en Ω, (1) donde: 1 L. J. ´ Alvarez-V´azquez, F. Javier Fern´andez, R. Mu˜noz-Sola l= 1, . . . , Ne, Ω⊂R3es un dominio con frontera ∂Ω suficientemente regular, Q= [0, T]×Ω y Σ = [0, T]×∂Ω, cl(t, x, u,∇ul) = PNe k,m=1 cl k,m(t, x, uk,∇um) + PNe k,m=1 cl 0,k,m(t, x, uk, um). Dentro del sistema de ecuaciones anterior consideraremos los dos subtipos siguientes: Simplificaci´on Cuasilineal. al i(t, x, ul,∇ul) = Pn j=1 al i,j(t, x)∂ul ∂xj+al i,0(t, x, ul),(2) cl k,m(t, x, uk,∇um) = δk,m Pn j=1 Ψl 1,0,m,j(t, x, ∂um ∂xj),(3) cl 0,k,m(t, x, uk, um) = Ψl 0,k,0,0(t, x, uk)+Ψl 0,k,m,0(t, x, uk)um.(4) Para esta clase de sistemas diremos que ul∈W1,2,2(I;H1(Ω), H1(Ω)∗)1∩C(I;L2(Ω)) con ul(0, x) = u0l(x) y l= 1, . . . , Ne, es una soluci´on d´ebil si: PNe l=1 ©<dul dt , vl>+RΩal(t, x, ul,∇ul)· ∇vl(x) + cl(t, x, u,∇u)vl(x)dx+ RΓl 1bl(t, x, ul)vl(x)dγª=PNe l=1 ©RΩgl(t, ·)vldx +RΓl 1hl(t, ·)vldγªc.p.t. t ∈I, ∀v∈H1(Ω)Ne. (5) Simplificaci´on Semilineal. al i(t, x, ul,∇ul) = Pn j=1 al i,j(t, x, ul)∂ul ∂xj+al i,0(t, x, ul),(6) cl k,m(t, x, uk,∇um) = Pn j=1 Ψl 1,k,m,j(t, x, uk)∂um ∂xj,(7) cl 0,k,m(t, x, uk, um) = Ψl 0,k,0,0(t, x, uk)+Ψl 0,k,m,0(t, x, uk)um.(8) Para esta clase de sistemas diremos que ul∈L2(I;H1(Ω)) ∩ Cw(I;L2(Ω)) con l= 1, . . . , Nles una soluci´on muy d´ebil cuando: PNe l=1 ©RΩul(T, x)vl(T, x)dx +RQal(t, x, ul,∇ul)· ∇vl+cl(t, x, u,∇u)vl−∂vl ∂t uldxdt+ RΣl 1bl(t, x, ul)vldγdtª=PNe l=1 ©RQglvldxdt +RΣl 1hlvldγdt +RΩu0l·vl(0,·)dxª, ∀v∈W1,∞,∞(I; [W1,∞(Ω)]Ne,([H1(Ω)]Ne)∗). (9) 1W1,p,q(I;V, V ∗) = {u∈Lp(I;V) : du dt ∈Lq(I;V∗)}. 2 Problemas de control en procesos de eutrofizaci´on La existencia de soluci´on de esta clase de sistemas se obtiene empleando t´ecnicas similares a las descritas en [3]. La idea b´asica para demostrar la existencia de soluci´on d´ebil consiste en imponer hip´otesis sobre los coeficientes para garantizar una serie de condiciones sobre el operador asociado al sistema tales como la acotaci´on y la pseudomonoton´ıa. Por otro lado, para garantizar la existencia de soluciones muy d´ebiles, ser´a necesario imponer restricciones sobre los coeficientes que nos permitan asegurar propiedades sobre el operador como la acotaci´on y la continuidad d´ebil en ciertos espacios. Los modelos que gobiernan los procesos de eutrofizaci´on est´an basados en ecuaciones en derivadas parciales de una gran complejidad debido a la variedad de fen´omenos que intervienen. En este trabajo hemos considerado un modelo relativamente simple en donde s´olo se consideran las siguientes cinco especies ([1]): 1. Fitoplancton (u1), 2. Nutriente Gen´erico, por ejemplo, el Nitr´ogeno (u2), 3. Zooplancton (u3), 4. Ox´ıgeno Disuelto (u4), 5. Detritos Org´anicos (u5). La interacci´on entre ellas en un dominio Ω ⊂R3(que vamos a suponer fijo) se puede modelar mediante el siguiente sistema de ecuaciones en derivadas parciales:                                                  ∂u1 ∂t +w·u1−∇·(µ1∇u1)−Au2 KN+u2u1+ (Kr+Kmf )u1+Kzu1 KF+u1u3=g1en Q, ∂u2 ∂t +w·u2−∇·(µ1∇u1) + CncAu2 KN+u2u1−CncKru1−CncKrdΘθ−20u5=g2en Q, ∂u3 ∂t +w·u3−∇·(µ1∇u1)−CfzKzu1 KF+u1u3+Kmzu3=g3en Q, ∂u4 ∂t +w·u4− ∇ · (µ1∇u1)−CocAu2 KN+u2u1+CocKru1+CocKrdΘθ−20u5=g4en Q, ∂u5 ∂t +w·u5− ∇ · (µ1∇u1)−Kmf u1−Kmzu3+KrdΘθ−20u5−Wfd ∂u5 ∂z =g5en Q, i= 1, . . . 5, µi∇ui·n+αiui=hisobre Σ, i= 1, . . . 5, ui(0) = u0ien Ω. (10) donde wes la velocidad del agua, Coc es el coeficiente estequiom´etrico ox´ıgeno-carbono, Cnc es el coeficiente estequiom´etrico nitr´ogeno-carbono, Krd es la tasa de regeneraci´on de los detritos, Θ es la constante de regeneraci´on t´ermica de los detritos, θes la temperatura del agua, Kres la tasa de respiraci´on end´ogena del fitoplancton, Kmf es la tasa de mortalidad del fitoplancton, Kmz es la tasa de mortalidad del zooplancton, Kzes el coeficiente de predaci´on del zooplancton, Wfd es la velocidad de sedimentaci´on de los detritos org´anicos, KFes la constante de saturaci´on media del fitoplancton, KNes la constante de saturaci´on 3 L. J. ´ Alvarez-V´azquez, F. Javier Fern´andez, R. Mu˜noz-Sola media del nitr´ogeno, µi,i= 1, . . . , 5 son los coeficientes de difusi´on de las especies y Aes la funci´on de luminosidad: A=µCθ−20 t I0 Is e−(φ1+φ2u2)z,(11) con I0la intensidad de luz incidente, Isla saturaci´on de luz, Ctla constante de crecimiento t´ermico del fitoplancton, φ1la atenuaci´on de la luz debida a la profundidad del agua, φ2 es la atenuaci´on de la luz debida a la presencia de fitoplancton y µla tasa m´axima de crecimiento para el fitoplancton. Si suponemos que φ2= 0 (el caso distinto de cero es an´alogo y s´olo requiere unas peque˜nas modificaciones en el sistema que hemos presentado al comienzo de este art´ıculo) el sistema anterior se puede considerar como un caso particular del sistema (1) y, dependiendo de la regularidad del t´ermino de convecci´on y los datos, lo incluiremos como un caso particular de la simplificaci´on cuasilineal o un caso particular de la simplificaci´on semilineal. En concreto, suponiendo que los t´erminos fuente y de frontera son nulos y que las condiciones iniciales tienen la regularidad requerida en cada caso, si: w·∇u∈L2(I;H1(Ω)∗), ∀u∈L2(I;H1(Ω)), se incluye en la simplificaci´on cuasilineal, w· ∇u∈L1+²∗(I;W1,q(Ω)∗), ∀u∈L2(I;H1(Ω)), con ²∗>0 y 1 < q < ∞, se incluye en la simplificaci´on semilineal. En este trabajo supondremos que la velocidad w, los t´erminos fuente gi, de contorno hiyαien (10) son cero y que las condiciones iniciales son suficientemente regulares (ui 0∈W7 5,10 3(Ω)). Bajo estas hip´otesis tiene sentido la formulaci´on d´ebil del problema y, adem´as, puede obtenerse regularidad adicional (ui∈W2,1 10 3 (Q)2∩C(Q) ) empleando t´ecnicas y resultados de [2]. Una vez demostrada la regularidad adicional para la soluci´on, podemos demostrar la unicidad de la misma, la continuidad y la derivabilidad con respecto a las condiciones iniciales. 2. El problema de control. Depuracion de aguas en Biorreactores. La idea fundamental de los biorreactores consiste en estancar el agua rica en nitr´ogeno en grandes tanques a los cuales a˜nadiremos una cierta cantidad de fitoplancton que dejaremos desarrollarse para que absorba el nitr´ogeno del agua. En el problema de control que planteamos tenemos dos tanques de igual volumen. El agua rica en nitr´ogeno llegar´ıa al primer tanque, al cual a˜nadir´ıamos una cierta cantidad de fitoplancton que dejar´ıamos crecer para que hiciese descender los niveles de nitr´ogeno hasta un cierto umbral m´aximo. En este primer tanque nos interesar´ıa tambi´en que se produjera una cierta cantidad de detritos org´anicos pues ser´ıan una buena fuente de abono agr´ıcola. Alcanzados los niveles m´aximo de nitr´ogeno y m´ınimo de detritos, filtramos los 2W2,1 p(Q) = {u∈Lp(Q) : ||u||Lp(Q)+||∂u ∂t ||Lp(Q)+P3 j=1 || ∂u ∂xi||Lp(Q)+P3 i,j=1 || ∂2u ∂xi∂xj||Lp(Q)<∞} 4 Problemas de control en procesos de eutrofizaci´on detritos y pasamos el agua a un segundo tanque para volver a repetir la operaci´on. El agua que salga de este segundo tanque ser´a agua pobre en nitr´ogeno, pero sin embargo, puede contener una gran cantidad de fitoplancton que tambi´en nos interesar´ıa minimizar. Por lo tanto, se trata de controlar la cantidad de fitoplancton que tenemos que a˜nadir a cada uno de los dos tanques, para que los niveles de nitr´ogeno est´en por debajo de unos m´aximos y los de detritos por encima de unos m´ınimos y de forma que la concentraci´on final de fitoplancton sea lo m´as reducida posible. Desde el punto de vista matem´atico este problema se puede plantear como un problema de control ´optimo con restricciones sobre el estado, en donde el control (ρ1, ρ2) ser´a la concentraci´on de fitoplancton que a˜nadiremos a cada uno de los tanques. Ecuaciones de estado. Tenemos dos sistemas de ecuaciones de estado, uno para cada tanque (Ω1y Ω2). En principio, vamos a prescindir de la ecuaci´on del ox´ıgeno disuelto pues no juega un papel relevante. Al cambiar de tanque, supondremos que el agua se mezcla y es por ello que las condiciones iniciales para el tanque 2 ser´an la media de las concentraciones, en el tiempo final, del tanque 1. •Tanque 1. Las variables de estado para el tanque 1 ser´an: u1,1(nitr´ogeno), u2,1(fitoplancton), u3,1(zooplancton) y u5,1(detritos org´anicos). El tiempo de permanencia en el primer tanque ser´a T1, siendo Q1= [0, T1]×Ω1y Σ1= [0, T1]×∂Ω1.                                                          ∂u1,1 ∂t − ∇ · (µ1∇u1,1)−Au2,1 KN+u2,1u1,1+ (Kr+Kmf )u1,1+Kzu1,1 KF+u1,1u3,1= 0 en Q1, ∂u2,1 ∂t − ∇ · (µ2∇u2,1) + CncAu2,1 KN+u2,1u1,1−CncKru1,1−CncKrdΘθ−20u5,1= 0 en Q1, ∂u3,1 ∂t − ∇ · (µ3∇u3,1)−CfzKzu1,1 KF+u1,1u3,1+Kmzu3,1= 0 en Q1, ∂u5,1 ∂t − ∇ · (µ5∇u5,1)−Kmf u1,1−Kmzu3,1+KrdΘθ−20u5,1−Wfd ∂u5,1 ∂z = 0 en Q1, i= 1,2,3,5,∇ui,1·n= 0 sobre Σ1, i= 1,3, ui,1(0) = ui,1 0en Ω1, u2,1(0) = u2,1 0+ρ1en Ω1, u5,1(0) = 0, en Ω1. (12) •Tanque 2. Las variables de estado para el tanque 1 ser´an: u1,2(nitr´ogeno), u2,2(fitoplancton), u3,2(zooplancton) y u5,2(detritos org´anicos). El tiempo de permanencia en el segundo tanque ser´a T2, siendo Q2= [0, T2]×Ω2y Σ2= [0, T2]×∂Ω2. 5 L. J. ´ Alvarez-V´azquez, F. Javier Fern´andez, R. Mu˜noz-Sola                                                          ∂u1,2 ∂t − ∇ · (µ1∇u1,2)−Au2,2 KN+u2,2u1,2+ (Kr+Kmf )u1,2+Kzu1,2 KF+u1,2u3,2= 0 en Q2, ∂u2,2 ∂t − ∇ · (µ2∇u2,2) + CncAu2,2 KN+u2,2u1,2−CncKru1,2−CncKrdΘθ−20u5,2= 0 en Q2, ∂u3,2 ∂t − ∇ · (µ3∇u3,2)−CfzKzu1,2 KF+u1,2u3,2+Kmzu3,2= 0 en Q2, ∂u5,2 ∂t − ∇ · (µ5∇u5,2)−Kmf u1,2−Kmzu3,2+KrdΘθ−20u5,2−Wfd ∂u5,2 ∂z = 0 en Q2, i= 1,2,3,5,∇ui,2·n= 0 sobre Σ2, i= 1,3, ui,2(0) = 1 µ(Ω1)RΩ1ui,1(T1)dx en Ω2, u2,2(0) = 1 µ(Ω1)RΩ1u2,1(T1)dx +ρ2en Ω2, u5,2(0) = 0 en Ω2. (13) Funcional de coste. Tal y como hemos explicado antes, nos interesa minimizar la concentraci´on de fitoplancton que sale del tanque 2: J(ρ1, ρ2) = 1 2||u2,2(T2)||2 L2(Ω2).(14) Restricciones sobre el estado. La concentraci´on final de nitr´ogeno en cada uno de los tanques debe de estar por debajo de unos niveles m´aximos y la concentraci´on final de detritos por encima de unos niveles m´ınimos:      B1(ρ1, ρ2) = 1 µ(Ω1)RΩ1u1,1(T1)dx ≤σ1, B2(ρ1, ρ2) = 1 µ(Ω2)RΩ2u1,2(T2)dx ≤σ2, B3(ρ1, ρ2) = 1 µ(Ω1)RΩ1u5,1(T1)dx ≥θ1, B4(ρ1, ρ2) = 1 µ(Ω2)RΩ2u5,2(T2)dx ≥θ2. (15) Por lo tanto, la formulaci´on del problema de control ´optimo ser´a la siguiente: (P) inf ©J(ρ1, ρ2)t.q. (ρ1, ρ2)∈ Uad, con (u1,1, u5,1, u1,2, u5,2)verificando (15)ª,(16) donde Uad ⊂W7 5,10 3(Ω1)×W7 5,10 3(Ω2) es un subconjunto acotado, convexo, cerrado y no vac´ıo. Podemos probar que este problema de control admite soluci´on no necesariamente ´unica. Adem´as, si (ρ1 0, ρ2 0) es una soluci´on del problema de control, entonces existir´a una constante γ≥0, −→ λ∈R4y elementos (u1,1 0, u2,1 0, u3,1 0, u5,1 0), (p1,1 0, p2,1 0, p3,1 0, p5,1 0)∈W2,1 10 3 (Q1)∩ C(Q1), (u1,2 0, u2,2 0, u3,2 0, u5,2 0), (p1,2 0, p2,2 0, p3,2 0, p5,2 0)∈W2,1 10 3 (Q2)∩ C(Q2) ´unicos tales que: (u1,1 0, u2,1 0, u3,1 0, u5,1 0)y(u1,2 0, u2,2 0, u3,2 0, u5,2 0) verifican las ecuaci´ones de estado (12)-(13), −→ λ∈∂κE(B(ρ1 0, ρ2 0)), donde E= [0, σ1]×[0, σ2]×[θ1,∞)×[θ2,∞), 6 Problemas de control en procesos de eutrofizaci´on (p1,1 0, p2,1 0, p3,1 0, p5,1 0) y (p1,1 0, p2,1 0, p3,1 0, p5,1 0) verifican las ecuaciones del estado adjunto:                                                                  ∂p1 0 ∂t −div(Aµ∇p1 0) = DE(c1 0)Tp1 0en Q1, ∂p2 0 ∂t −div(Aµ∇p2 0) = DE(c2 0)p2 0en Q2, ∇p1 0·n= 0 sobre Σ1, ∇p2 0·n= 0 sobre Σ2, p1,1 0(T1) = λ1+1 µ(Ω1)RΩ2p1,2 0(0)dx en Ω1, p2,1 0(T1) = 1 µ(Ω1)RΩ2p2,2 0(0)dx en Ω1, p3,1 0(T1) = 1 µ(Ω1)RΩ2p3,2 0(0)dx en Ω1, p5,1 0(T1) = λ3en Ω1, p1,2 0(T2) = λ2, p5,2 0(T2) = λ4, pi,2 0(T2) = 0, i = 2,3, en Ω2, (17) donde ui= (u1,i, u2,i, u3,i, u5,i) y DE(ui) es la siguiente matriz:            −CncAKN (KN+u1,i)2u2,i −CncAu1,i KN+u1,i +CncKr0CncKrdΘθ−20 AKN (KN+u1,i)2u2,i Au1,i KN+u1,i −(Kr+Kmf )−KzKF (KF+u2,i)2u3,i −Kzu2,i KF+u2,i 0 0CfzKzKF (KF+u2,i)2u3,i CfzKzu2,i KF+u2,i −Kmz 0 0Kmf Kmz −KrdΘθ−20            (18) Se satisface la siguiente condici´on de optimalidad: γJ0(ρ1 0, ρ2 0)(ρ1−ρ1 0, ρ2−ρ2 0)+ 1 µ(Ω1)RΩ1(ρ1−ρ1 0)p2,1dx +1 µ(Ω2)RΩ2(ρ2−ρ2 0)p2,2dx ≥0,∀(ρ1, ρ2)∈ Uad. (19) 3. Resoluci´on num´erica del problema. Un algoritmo de puntos interiores. La resoluci´on num´erica del problema de control que hemos planteado, presenta una serie de complejidades como son la resoluci´on de la ecuaci´on de estado, la resoluci´on de las ecuaciones linealizadas (necesarias para el c´alculo del gradiente del funcional de coste 7 L. J. ´ Alvarez-V´azquez, F. Javier Fern´andez, R. Mu˜noz-Sola y de las restricciones) o la propia resoluci´on del problema de minimizaci´on asociado. Para la resoluci´on num´erica de la ecuaci´on de estado hemos empleado una discretizaci´on por elementos finitos combinada con un algoritmo de punto fijo para resolver la no linealidad y, en el caso de las ecuaciones linealizadas, hemos empleado un algoritmo de relajaci´on para tratar el acoplamiento entre especies. Por ´ultimo, la resoluci´on del problema de minimizaci´on asociado, la hemos realizado empleando el c´odigo IPOPT, un algoritmo de puntos interiores cuya descripci´on se puede ver en [4]. A continuaci´on presentamos, a modo de ejemplo, la evoluci´on del algoritmo en un test con par´ametros acad´emicos en donde se alcanza saturaci´on en una de las restricciones: iter objective inf_pr inf_du lg(mu) ||d|| 0 3.6024684e+01 0.00e+00 5.09e+00 0.0 0.00e+00 1 1.8439766e+01 3.10e-01 1.74e+00 0.2 4.88e+00 2 1.8474918e+01 9.91e-04 5.75e-01 -1.1 3.45e-01 3 1.8057410e+01 1.65e-03 9.65e-02 -2.3 1.94e-01 4 1.8055934e+01 3.79e-04 3.56e-01 -3.2 2.01e-03 5 1.8055104e+01 2.21e-07 2.37e-04 -4.7 1.62e-03 6 1.8055061e+01 7.15e-11 1.13e-06 -9.0 6.92e-06 7 1.8055061e+01 2.66e-13 4.33e-12 -9.0 1.52e-09 En las siguientes figuras se puede apreciar el comportamiento de las concentraciones de nitr´ogeno en el tiempo final para cada uno de los tanques: Tanque 1 Tanque 2 Referencias [1] M. Drago, B. Cescon, L. Iovenitti, A three-dimensional numerical model for eutrophication and pollutant transport. Ecological Modelling, 145 (2001), 17-34. [2] O. A. Ladyzenskaja, V. Solonnikov, N. N. Ural´ceva, Linear and Quasilinear Equations of Parabolic Type. Translations of Mathematical Monographs, Vol. 23. A.M.S., Providence, 1968. [3] T. Roub´ıcek, Nonlinear Partial Differential Equations with Applications. Birkh¨auser, Berlin, 2005. [4] A. W¨achter, L. T. Biegler, On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. Math. Program., Ser. A, 106 (2006), 25-27. 8