Modelado, simulación y control del uso del agua en la agricultura

Modelado, simulación y control del uso del agua en la agricultura Javier Arriaga Sevilla Abril 2015 UNIVERSIDAD DE SEVILLA ESCUELA TÉCNICA SUPERIOR DE INGENIERÍA TESIS DOCTORAL Modelado, simulación y control del uso del agua en la agricultura Tesis presentada ante el Departamento de Ingeniería de Sistemas y Automática, Escuela Técnica Superior de Ingeniería, como parte de los requisitos para la obtención del título de Doctor por la Universidad de Sevilla. de Javier Arriaga Sevilla Sevilla, Abril 2015 Directores: Francisco Rodríguez Rubio José Luis Muriel Fernández Ivan Francisco García Tejero A mis padres Agradecimientos La elaboración de una tesis doctoral es, desde sus inicios, una labor ardua en la que confluyen multitud de piezas en movimiento; hubiera sido imposible llevarla a cabo sin la ayuda de muchas personas, ayuda de la que deseo dar fe en las próximas líneas. En primer lugar, me gustaría dar las gracias a mis padres. Sin su amor, sin su ayuda y sin su santa paciencia, hubiera sido imposible llegar hasta este punto. Lo mismo se puede decir del resto de mi familia, que siempre me animó y espoleó para lograr mi objetivo. Debo expresar la mayor gratitud a mis tres directores de tesis. A Francisco R. Rubio, por no haber tirado la toalla conmigo a pesar de que en muchas ocasiones le di motivos más que sobrados para ello; a José Luis Muriel, por haberme animado de forma constante a seguir adelante y dar estructura al trabajo; y por supuesto, a mi amigo Iván García porque en alguna ocasión casi literalmente me levantó cuando yo mismo estaba dispuesto a abandonar. No puedo olvidarme de mis compañeros de departamento en el centro IFAPA de Alcalá del Río: Karl, con el que di los primeros y titubeantes pasos de esta tesis; Mª Antonia y Natalia, siempre presentes y a veces imprescindibles en mi día a día; Aura, Antonio, Victor Hugo, Pepe, Gonzalo, Manolo y tantos otros. Ni de los muchos amigos que he hecho en el centro, amigos que han hecho estos años una experiencia maravillosa: Jorge, Miguel Ángel, Juan Antonio, Almudena, Manuel, Ana María, las hermanas Gómez, Carmen y las maravillosas Nieves y Ana. En estos años he trabajado en varios lugares, pero si ha habido una constante ha sido la academia, donde Manuel ha sido más que un jefe, un amigo. Gracias. Finalmente, no puedo dejar de mencionar a mi mejor amiga, una persona a la que conocí a los dieciocho años y que se ha convertido en una parte imprescindible de mi vida. Gracias, María Dolores, porque siempre confiaste en v vi Agradecimientos mí y porque siempre pude contar contigo cuando lo necesité; y te he necesitado mucho más a menudo de lo que me gusta reconocer. Por supuesto, es necesario agradecer al Ministerio de Educación y Ciencia y al Instituto Nacional de Innovación Agraria (INIA) por la beca FPI con la que se financió el desarrollo del trabajo de mi tesis. Además de a la Junta de Andalucía por el proyecto P09-AGR-4782 «El Papel de la Humedad del Suelo en la Conservación de Suelo, Agua y Carbono en Cuencas Agrícolas» y al Ministerio de Ciencia e Innovación por el proyecto del Plan Nacional de I+D AGL200912936-C03-02 «Integración de medidas en cuencas agrícolas a diferentes escalas en un modelo de agua y sedimentos.» Javier Arriaga Sevilla Sevilla, Abril de 2015. Índice general Agradecimientos v Índice de figuras xi Índice de tablas xv Resumen xvii Summary xix 1. Introducción 1 1.1.Antecedentes.............................. 1 1.2.Objetivos ................................ 4 1.3. Estructura de la Tesis . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.4. Publicaciones relacionadas . . . . . . . . . . . . . . . . . . . . . . 8 1.4.1. Artículos ............................ 8 1.4.2. Presentaciones en congresos . . . . . . . . . . . . . . . . . 9 Parte I Modelado del agua en suelo 11 2. Calibración por Redes Neuronales Artificiales 13 2.1.Introducción.............................. 15 2.1.1. Sensores electromagnéticos para la medición de agua en suelo .............................. 15 2.1.2. Redes Neuronales Artificiales . . . . . . . . . . . . . . . . 19 2.2.Materialymétodos .......................... 21 2.2.1. Ámbito del ensayo . . . . . . . . . . . . . . . . . . . . . . . 21 2.2.2. Recogida y tratamiento de muestras . . . . . . . . . . . . . 22 2.2.3. Calibración........................... 26 2.3. Resultados y discusión . . . . . . . . . . . . . . . . . . . . . . . . 28 2.3.1. Recopilación de datos . . . . . . . . . . . . . . . . . . . . . 28 vii viii Índice general 2.3.2. Cálculo de las redes neuronales artificiales. . . . . . . . . 30 2.3.3. Evaluación de resultados . . . . . . . . . . . . . . . . . . . 33 2.3.4. Modelo de tres sensores. . . . . . . . . . . . . . . . . . . . 37 2.4.Conclusiones.............................. 40 3. Modelo del contenido de humedad en suelo 41 3.1.Introducción.............................. 41 3.2. Material y métodos . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.2.1. El modelo del cubo . . . . . . . . . . . . . . . . . . . . . . 44 3.2.2. Modificaciones......................... 47 3.2.3. Entorno de simulación . . . . . . . . . . . . . . . . . . . . 51 3.2.4. Calibración........................... 52 3.3.Resultados ............................... 53 3.4.Conclusiones.............................. 57 4. Control del contenido de agua en suelo 59 4.1.Introducción.............................. 59 4.2. Material y métodos . . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.2.1. Determinación del punto de operación . . . . . . . . . . . 62 4.2.2. Diseño de controladores . . . . . . . . . . . . . . . . . . . 63 4.3.Resultados ............................... 67 4.3.1. Linealización del sistema . . . . . . . . . . . . . . . . . . . 67 4.3.2. Control Proporcional . . . . . . . . . . . . . . . . . . . . . 70 4.3.3. Control Proporcional Integral (PI) . . . . . . . . . . . . . . 70 4.3.4. Control Proporcional Integral Derivativo (PID) . . . . . . 72 4.3.5. Simulación de controladores . . . . . . . . . . . . . . . . . 73 4.4.Conclusiones.............................. 77 Parte II Modelado de los efectos del Riego Deficitario sobre cultivos leñosos 79 Introducción................................. 82 Objetivos................................... 86 5. Efectos del estrés hídrico a largo plazo 87 5.1. Material y métodos . . . . . . . . . . . . . . . . . . . . . . . . . . 87 5.1.1. Parcela experimental . . . . . . . . . . . . . . . . . . . . . 87 5.1.2. Tratamientos de RD . . . . . . . . . . . . . . . . . . . . . . 88 5.1.3. Análisis estadístico y numérico; modelo analítico . . . . . 89 5.2. Resultados y discusión . . . . . . . . . . . . . . . . . . . . . . . . 95 5.2.1. Relaciones hídricas y respuesta fisiológica del cultivo . . . 95 Modelado, simulación y control del uso del agua en la agricultura ix 5.2.2. Respuesta productiva del cultivo a riego deficitario sostenido y funciones de productividad del agua anuales. . . . 96 5.2.3. Funciones de pérdida de producción . . . . . . . . . . . . 99 5.3.Conclusiones.............................. 105 6. Efectos del estrés hídrico a medio plazo 107 6.1.Materialymétodos .......................... 107 6.1.1. Parcela experimental . . . . . . . . . . . . . . . . . . . . . 107 6.1.2. Diseño experimental. Tratamientos RDC . . . . . . . . . . 108 6.1.3. Modelando la función de producción del cultivo . . . . . 110 6.1.4. Modeloaditivo......................... 111 6.1.5. Modelo por ajuste con términos cruzados . . . . . . . . . . 112 6.2. Resultados y discusión . . . . . . . . . . . . . . . . . . . . . . . . 116 6.2.1. Producción y tratamientos de RDC . . . . . . . . . . . . . 116 6.2.2. Ajuste sin términos cruzados . . . . . . . . . . . . . . . . . 118 6.2.3. Ajuste con términos cruzados . . . . . . . . . . . . . . . . 120 6.3.Conclusiones.............................. 126 7. Planificación de estrategias de riego 129 7.1.Introducción.............................. 129 7.2.Materialymétodos .......................... 132 7.2.1. Planificación a largo plazo . . . . . . . . . . . . . . . . . . 132 7.2.2. Planificación a medio plazo . . . . . . . . . . . . . . . . . 140 7.3.Resultados ............................... 141 7.3.1. Planificación a largo plazo . . . . . . . . . . . . . . . . . . 141 7.3.2. Planificación a medio plazo . . . . . . . . . . . . . . . . . 152 7.4. Propuesta de aplicación mixta . . . . . . . . . . . . . . . . . . . . 154 7.4.1. Hipótesis de partida . . . . . . . . . . . . . . . . . . . . . . 154 7.4.2. Metodología .......................... 155 7.5.Conclusiones.............................. 157 8. Riego en almendros 159 8.1.Introducción.............................. 160 8.2.Materialymétodos .......................... 161 8.2.1. Condiciones del experimento . . . . . . . . . . . . . . . . 161 8.2.2. Tratamientos de riego . . . . . . . . . . . . . . . . . . . . . 163 8.2.3. Mediciones en planta . . . . . . . . . . . . . . . . . . . . . 163 8.2.4. Análisis de datos y modelado . . . . . . . . . . . . . . . . 164 8.2.5. Diseño del controlador . . . . . . . . . . . . . . . . . . . . 165 8.3. Resultados y discusión . . . . . . . . . . . . . . . . . . . . . . . . 167 8.3.1. Modelo multivariable LTI en tiempo discreto . . . . . . . 168 xvi Índice de tablas 6.4. Tasas de pérdida de producción (kg/kg) medidas (med.) y previstas (prev.) con respecto a la producción del tratamiento control según tratamiento y año. . . . . . . . . . . . . . . . . . . . . . . . 125 7.1. Errores del modelo según tratamiento . . . . . . . . . . . . . . . . 142 7.2. Distribución de riegos óptima (en mm): resultados de aplicar año a año la secuencia calculada para RDS con RDC. Las pérdidas se miden en tanto por uno (kg/kg). El reparto óptimo es: 79 mm en floración, 187 mm en crecimiento y 51 mm en maduración. . . . 156 Resumen La agricultura en Europa, y más concretamente en regiones periféricas como Andalucía, supone un elemento fundamental para el desarrollo económico, especialmente en las zonas rurales, actuando como elemento vertebrador y de cohesión tanto a nivel social como productivo. Dentro de las innumerables limitaciones con las que se encuentra el agricultor, la disponibilidad de agua es quizá una de las más determinantes, no sólo desde el punto de vista de la disponibilidad, sino también de la incertidumbre que genera un contexto de cambio climático global con una creciente demanda que amenaza con sobrepasar los recursos disponibles. Esta incertidumbre resulta patente tanto en los sistemas productivos de secano como en los de regadío, presentándose como necesaria la introducción de diferentes herramientas y modelos que permitan alcanzar una gestión más eficiente y sostenible de los recursos hídricos, además de la posibilidad de plantear distintos escenarios futuros y elaborar respuestas adecuadas adelantándose a los acontecimientos. El presente trabajo de tesis doctoral ofrece diferentes enfoques de la problemática de la escasez de agua, en base a: i) la introducción de diferentes sensores que posibiliten el modelado y control del sistema Agua-Suelo-Planta-Atmósfera al ofrecer medidas de distintas variables de los componentes del sistema. ii) el planteamiento de modelos matemáticos y computacionales que permitan una caracterización y representación de una realidad compleja como es el citado sistema. iii) la propuesta de varios sistemas de planificación y control automático de distintas variables de interés (el contenido de agua, su productividad, la producción del cultivo y otras variables relacionadas con la fisiología del cultivo). Para ello, los trabajos experimentales se han desarrollado en un sistema de agricultura extensiva de secano, representativa de casi un 70% de los agroecosistemas andaluces, y en dos cultivos leñosos de regadío que cuentan con una elevada presencia en nuestra Comunidad Autónoma: el cultivo de cítricos y el del almendro. A partir de los resultados obtenidos se puede concluir que xvii xviii Resumen es posible encontrar representaciones matemáticas relativamente sencillas de procesos naturales muy complejos y aparentemente caóticos. Estas representaciones han permitido en gran medida modelar el comportamiento del agua en el suelo en condiciones de secano, y comprobar que dichos modelos eran capaces de recoger las características inherentes de cada sistema, como por ejemplo el sistema de laboreo empleado, la profundidad del horizonte radicular o las propiedades físicas del medio. Por otra parte, en cuanto a los sistemas de regadío, los modelos obtenidos permiten extraer conclusiones no sólo sobre la influencia del aporte de agua en la producción final de los cultivos, sino también sobre la capacidad de adaptación de dichos cultivos a situaciones prolongadas de escasez. Estos modelos permiten a su vez discriminar el potencial efecto de un determinado nivel de estrés hídrico sobre los diferentes periodos fenológicos del cultivo; además, permiten concluir que el tiempo necesario para la adaptación de un cultivo y la estabilización de sus producciones a una situación prolongada de estrés hídrico depende tanto del tiempo como del nivel de estrés aplicado a lo largo de sucesivas campañas. Asimismo, los modelos hacen posible la planificación a largo plazo de estrategias de riego deficitario y sistemas de control automático para optimizar la implementación de dichas estrategias y maximizar los valores de producción final. En este sentido, se propone el empleo de técnicas clásicas de control como los controladores Proporcional-Integral-Diferencial o control mediante observador por filtro de Kalman. Teniendo en cuenta el conjunto de resultados obtenidos en este trabajo, se puede concluir que el empleo de técnicas de ingeniería de modelado, simulación y control resulta una alternativa viable para el mundo agrícola, y que sin duda puede ofrecer respuestas a las diferentes problemáticas del sector. Summary Agriculture in Europe, and especially in peripheral regions like Andalusia, is a foundation for economic development, especially in rural areas, acting as the backbone and cohesion element for society and economy. Among the many constraints the farmer finds himself under, water availability is one of the most, not just in regard to availability, but also to the uncertainty generated by a context of global climate change and a growing demand that threatens to overcome the available resources. This uncertainty is apparent both in rain fed and irrigated production systems, stating the necessity of introducing different tools and models to achieve a more efficient and sustainable management of water resources, raising the possibility of different future scenarios and the development of appropriate responses beforehand. The present Ph.D. Thesis offers different approaches to the problem of water scarcity, based on: i) the introduction of different sensors that makes possible the modelling and control of the Water-Soil-Plant-Atmosphere system through providing information about several variables and system components. ii) the approach to mathematical and computational models that allow the characterization and representation of complex realities as the abovementioned model. iii) the proposal of different planning and automatic control systems for several variables of interest (water content, productivity, crop production and other variables related to the crop physiology). In order to fulfill this objective, experimental work was developed in an extensive rain fed agricultural system, which represents almost a 70% of the andalusian agro-ecosystems and in two woody irrigated crops with a high presence in our region: citrus and almond trees. From the results obtained it can be concluded that it is possible to find simple mathematical representations of complex and seemingly chaotic natural processes. These representations allowed to model the behaviour of soil water under rain fed conditions and to attest that those models were able to take into account the system inherent characteristics, such as the tillage system used, the depth of the root horizon or the xix xx Summary soil physical properties. Moreover, in terms of irrigation systems, the obtained models allow to reach conclusions not just regarding the influence of water intake in the crop final production but also the crop adaptability to prolonged shortages. These models make possible to discriminate the potential effect of a given water stress level during different phenological stages; they also support the conclusion that the time needed by the crop to adapt itself and stabilize its production under prolonged water stress situations depends both on the time and the applied level of stress in the successive campaigns. The models also enable long-term planning of deficit irrigation strategies and automatic control systems to optimize the implementation of these strategies, thus maximizing the final productive value. In this regard, the use of conventional control techniques such as Proportional-Integral-Differential regulators and control using a Kalman filter observer are proposed. From the obtained results, it can be concluded that the use of engineering techniques as modelling, simulation and control is a viable option for the agricultural world and it can certainly provide answers to several problems of the sector. Capítulo 1 Introducción 1.1. Antecedentes El agua no es solamente el recurso natural más limitado que existe: es el combustible de la economía, motor de desarrollo y fuente de vida, tal vez incluso sinónimo de ésta. Ha sido el agua y su disponibilidad la que ha marcado históricamente la aparición, consolidación y decadencia de las primeras civilizaciones. La civilización egipcia floreció gracias al cultivo de los terrenos inundables del Nilo, y organizaba su vida pública en función de los ciclos de inundación de éste. Las civilizaciones del sudeste asiático, con China a la cabeza, ya concentraban a gran parte de la población mundial gracias al aprovechamiento del cultivo de regadío por inundación en la antigüedad. Incluso las civilizaciones consideradas como las más antiguas del mundo, la Sumeria y la Babilonia, surgieron en el célebre «creciente fértil» que se encontraba entre el río Tigris y el Eúfrates. Fue la sobreexplotación de los recursos hídricos en ese fértil territorio de la antigüedad, precisamente, lo que llevó a la desertización del territorio y a la inevitable decadencia de los imperios que en él se asentaban. Desde entonces, naciones e imperios han prosperado y han caído por múltiples motivos. En su ensayo «Colapso», Diamond [2005] describe cinco factores que llevaron al colapso a sociedades pasadas: 1 2 Capítulo 1. Introducción Daño ambiental. Cambio climático. Vecinos hostiles. Pérdida de comercio. Respuesta de la sociedad a problemas ambientales. De estos cinco factores, tres están fuertemente correlacionados con la disponibilidad de recursos hídricos. Obviamente el cambio climático y el daño ambiental afectan a otras variables (v.g., el ejemplo de la Isla de Pascua, donde la madera era la base de toda la tecnología de su civilización y sus habitantes cortaron hasta el último árbol). Sin embargo, prácticamente todas esas variables están a su vez fuertemente correlacionadas con el agua: la ausencia de árboles afecta a la capacidad del suelo para retener el agua, lo que provoca fuertes escorrentías en épocas de lluvias que arrastran tierra fértil consigo, produciendo desertización. El cambio climático, que históricamente siempre ha sido natural e irreversible (y esto último parece que lo sigue siendo), afecta inevitablemente a la disponibilidad de agua. Y la respuesta a los problemas ambientales acaba siendo lo que determina si las sociedades sobreviven, para lo cual han de responder a una pregunta: cómo gestionar el agua disponible. Según el informe de 2010 de las Naciones Unidas (DESA [2013]), revisado en 2012, se estima que la población humana en el año 1800 era de unos mil millones de personas. Doscientos años después, a pesar de haber vivido varias guerras mundiales, la población se había multiplicado por seis hasta los seis mil millones. Diez años después, había crecido hasta los siete mil millones. Este crecimiento exponencial, que se hace más acusado a partir de la década de los sesenta del pasado siglo, implica una necesidad de cultivo de alimentos que también ha de ser exponencialmente creciente, lo que implica un mayor uso de agua; o, al menos, un mejor uso de ella. Las proyecciones del citado informe indican tres escenarios posibles: una estabilización de la población en torno a los diez mil millones; una caída de la población que volvería a los seis mil doscientos millones a finales de este siglo; y la continuación de este crecimiento exponencial que nos llevaría a los dieciséis mil millones a final de siglo. En cualquiera de los tres escenarios, la presencia de un cambio climático cuyos efectos podrán ser mitigados pero en ningún caso neutralizados (Field y otros [2014]) afectará inequívocamente a la actividad agrícola y a la disponibilidad de agua. La figura 1.1 muestra esta evolución y tres proyecciones diferentes que muestran posibles futuros. Modelado, simulación y control del uso del agua en la agricultura 3 Figura 1.1: Proyecciones de población mundial según la ONU. Fuente: Loren Cobb para Wikimedia Project. Incluso en el mejor de los escenarios posibles, con la población estancada o en un lento retroceso, y con un cambio climático atemperado, el problema de la distribución de agua no es futuro: es presente. Según la FAO (Alexandratos y otros [2012]) la escasez de agua afecta a más del 40% de la población mundial, y las proyecciones de demanda de agua y alimentos para el futuro indican que el problema de la escasez no puede hacer sino aumentar. Por lo tanto, se puede concluir que el agua es el mayor tesoro de una civilización, y su gestión ha de ser una de las principales prioridades de cualquier gobierno. Con la población mundial experimentando un crecimiento exponencial es fundamental saber qué cantidad de agua se debe emplear en la producción 4 Capítulo 1. Introducción de cultivos agrícolas y, más importante aún si cabe, saber la cantidad de agua de la que se dispone. 1.2. Objetivos En concordancia con los antecedentes expuestos se planteó en este trabajo de Tesis como objetivo general el análisis y, en su caso, el modelado de los procesos implicados en el ciclo del agua a través del sistema continuo SueloPlanta-Atmósfera (SPA continuum). Ello permitiría una mejor comprensión de las relaciones entre estos componentes del sistema SPA a través de los cuales circula el agua, desde su almacenamiento en el suelo (procedente de la precipitación y/o riego), su absorción por la planta y su devolución a la atmósfera en forma de vapor por evaporación directa desde el suelo y a través de la transpiración de la planta. El conocimiento y modelado de estos procesos y la definición de sus relaciones contribuirán, sin duda, al diseño de estrategias sostenibles para la optimización del uso de los recursos hídricos, tanto en los sistemas productivos agrarios de secano bajo distintos tratamientos de manejo del suelo, como especialmente en los de regadío. Se pretende, como producto final, que el usuario disponga de nuevas tecnologías que racionarían de forma efectiva el uso del agua disponible y le faciliten la toma de decisiones a la hora de establecer sistemas de manejo del suelo y calendarios de riego que adapten la oferta real del recurso a la demanda exigida por los cultivos. Como primer objetivo parcial, se pretende la evaluación y modelado del estado hídrico en la zona no saturada del suelo en sus componentes espacial y temporal mediante el uso de Redes Neuronales Artificiales (RNA) como herramienta en ensayos de calibración de sensores. Una vez conseguida una calibración fiable de los sensores de humedad, pueden emplearse como fuentes de señales de información para su aplicación a modelos. La aplicación del modelo del cubo 1permitirá conocer la dinámica de la relación suelo-agua, y con ella el grado de disponibilidad para la planta, mediante el análisis de los procesos estocásticos implicados. Diversas modificaciones sobre este modelo permitirán su transformación en un modelo de parámetros distribuidos para aplicarlo a nivel de parcela; una vez conseguido este objetivo, se 1Introducido por Rodriguez-Iturbe [2000]. Modelado, simulación y control del uso del agua en la agricultura 5 pretende introducir en el mundo de la agronomía diversos sistemas de control clásico con el objetivo de estudiar su rendimiento en ensayos de distintos tipos de manejo. Una vez que el agua es captada desde el suelo por la planta, segundo componente vivo del sistema, se propone el estudio para establecer las necesidades hídricas reales de cultivos frutales de interés en Andalucía como cítricos y almendros, escogidos además por sus distintas características fenológicas y su comportamiento en cuanto a relaciones hídricas se refiere. Así, en los cítricos, el objetivo fue integrar en modelos analíticos la respuesta productiva a distintos tratamientos de riego deficitario, estableciendo las funciones de productividad del agua y las de pérdida de producción a largo y medio plazo. Mientras que en el caso del almendro se planteó realizar un análisis multifactorial del efecto de aportes diferenciales de agua sobre distintas variables fisiológicas mediante la aplicación de sistemas multivariables Lineales e Invariantes en el Tiempo (LTI) que permitieran identificar el estado hídrico óptimo de la planta en función de los aportes de riego recibidos. En resumen, esta Tesis plantea de forma secuencial los siguientes hitos: 1. Modelado del agua en el suelo Definir a partir de una red de sensores de humedad la evolución espacio-temporal del estado hídrico del suelo en función de los aportes naturales del agua por precipitación en distintos tratamientos de manejo del suelo: laboreo convencional y siembra directa. • Se incluye una propuesta de calibración de los sensores empleando Redes Neuronales Artificiales. Desarrollar un modelo dinámico observable y controlable que permita predecir con precisión el estado del agua en la zona no saturada del suelo y con ello la disponibilidad de la misma para el cultivo que soporta. • Se incluye el diseño y simulación de diversos sistemas de control clásico con el objetivo de estudiar su respuesta en ensayos de calibración de sensores. 2. Modelado de los efectos del Riego Deficitario sobre cultivos leñosos En cítricos Capítulo 2 Evaluación y calibración de sensores de humedad en campo y laboratorio usando redes neuronales. Resumen El estudio de la variabilidad espacial de la humedad del suelo a escala de parcela o cuenca agrícola requiere el uso de redes de sensores de humedad de bajo coste, que suelen mostrar una fiabilidad limitada y requieren de calibraciones específicas, especialmente en suelos con elevados contenidos en arcilla. Inmerso en el contexto de un estudio de la distribución espacial y temporal de la humedad del suelo en un ensayo a largo plazo de manejo de suelo a escala de parcela, el presente trabajo pretende plantear una calibración más fiable y eficiente de sondas de humedad mediante un análisis mixto campo-laboratorio. Para la calibración de campo se dispone de datos gravimétricos; para la calibración en laboratorio se emplearon columnas de suelo inalterado que tras ser saturadas fueron desecadas en un entorno controlado mediante una estufa de aire forzado, mientras se monitorizaba la evolución de su peso con básculas de 13 14 Capítulo 2. Calibración por Redes Neuronales Artificiales precisión y la de su humedad volumétrica medida con distintas sondas capacitivas de humedad Decagon. Tras obtener curvas de secado y la densidad aparente (gr/cm3)que relaciona las humedades gravimétrica y volumétrica, es posible realizar una calibración mejorada y específica para cada tipo de suelo. Las redes neuronales artificiales son particularmente útiles para el modelado de procesos físicos y el ajuste de dichos modelos; en este capítulo se propone el empleo de estas herramientas para obtener calibraciones para las sondas analizadas en el tipo de suelo objeto de estudio. Los resultados muestran que dichas calibraciones permiten mejorar la precisión de las mediciones de humedad realizadas. Una versión de este trabajo se presentó en las Jornadas de Estudio de la Zona No Saturada del Suelo de Salamanca, 2011: J. Arriaga, K. Vanderlinden, J.L. Muriel, M. Morón, M.G. Ortega, F.R. Rubio. «Calibración de sensores de humedad capacitivos usando redes neuronales.» Actas de las X Jornadas de Investigación de la Zona no Saturada del Suelo, Salamanca (España), 2011 I.S.B.N.: 978-84-694-6642-1 Modelado, simulación y control del uso del agua en la agricultura 15 2.1. Introducción La medición adecuada del contenido de agua en suelo resulta crítica para la estimación de los balances de agua y energía, así como para comprender los procesos biológicos y químicos en todo el sistema suelo-planta (Vereecken y otros [2008], Robinson y otros [2003]). Puede encontrarse una revisión de técnicas de medición en Hopmans y Simunek [1999], incluyendo técnicas de medida geofísicas y remotas. 2.1.1. Sensores electromagnéticos para la medición de agua en suelo En concreto, los dispositivos basados en el análisis de las reflexiones de ondas electromagnéticas se han abierto un considerable hueco en el mercado. Concretamente, el estudio de la reflectometría en el dominio del tiempo para la medida de agua en suelo a través de las propiedades eléctricas del medio fue propuesto por Topp y otros [1980]; desde entonces, la tecnología de medición ha ido mejorando hasta convertirse en el método más popular para la medida de humedad en suelo frente a técnicas como la penetración por radar, inducción electromagnética, sensores de pulso de calor y otros. Simultáneamente, las técnicas de medida basadas en sondas de capacitancia y reflectometría en el dominio de la frecuencia (Frequency Domain Reflectometry, FDR) se han convertido en una alternativa viable debido a sus inherentes ventajas (Kizito y otros [2008]): bajo coste, facilidad para almacenamiento de datos, posibilidad de monitorización continua y funcionamiento en un amplio rango de tipos de suelos. El fundamento teórico de estos dispositivos se basa en la teoría electromagnética clásica de Maxwell (Pozar [2009]). Un condensador es un dispositivo eléctrico formado por dos placas metálicas con un dieléctrico separándolas, tal y como muestra la figura 2.1. Aunque éste es el modelo mostrado habitualmente, no es el único posible: los condensadores empleados habitualmente para electrónica, por ejemplo, tienen forma cilíndrica, con una placa cilíndrica exterior y una interior. 16 Capítulo 2. Calibración por Redes Neuronales Artificiales Figura 2.1: Esquema básico de un condensador de placas planoparalelas. El dieléctrico de separación (puede ser incluso el aire o el vacío) tendrá unas ciertas propiedades dieléctricas: permitividad eléctrica ϵ, permeabilidad magnética µy conductividad σ. Considerando que la ecuación de carga y descarga del condensador permite obtener la capacidad C y que ésta puede obtenerse a partir de la geometría y las propiedades físicas del dieléctrico, la relación entre estas propiedades y la humedad en suelo puede establecerse a partir de Topp y otros [1980]. Para comprender el funcionamiento de estos sensores es necesario hacer una breve introducción a la reflectometría, esto es, el estudio de las reflexiones; para explicar este proceso, hay que remitirse al estudio de la propagación de ondas en una línea de transmisión. El diagrama básico de una línea de transmisión puede observarse en la figura 2.2. Figura 2.2: Diagrama básico de una línea de transmisión. Modelado, simulación y control del uso del agua en la agricultura 17 Se dice que una línea de transmisión está adaptada cuando al aplicar una señal (en forma de diferencia de potencial en la puerta de entrada), ésta no sufre reflexiones. Sin embargo, cuando la señal de información se ve reflejada, se dice que la línea no está adaptada, y esto se debe a que la impedancia de la carga (impedancia vista en el extremo final de la línea) no coincide con la impedancia característica de la línea. Por lo tanto, a partir de una impedancia de carga conocida, aplicando una señal eléctrica a la entrada de una línea de transmisión, es posible a partir del estudio de las reflexiones conocer la impedancia característica de la línea. Naturalmente, dicha impedancia es una propiedad que depende de las propiedades electromagnéticas del dieléctrico que aísla ambos hilos de la línea de transmisión (permitividad eléctrica y permeabilidad magnética), además de la geometría y frecuencia de operación. Por lo tanto, controlando la geometría y la frecuencia, es posible determinar las propiedades del dieléctrico y a partir de ellas determinar el contenido de agua en suelo. Dentro de este estudio de la reflectometría encontramos dos técnicas: 1. La reflectometría en el dominio de la frecuencia (FDR) es la técnica que usan las llamadas sondas de capacitancia. Se basan la emisión de señales de información a cierta frecuencia y, a partir del estudio de los resultados en régimen permanente determinar la impedancia. 2. La reflectometría en el dominio del tiempo (TDR), por el contrario, trabaja en el dominio del tiempo, por lo que se centra en el estudio de los transitorios. Se analiza el desfase temporal entre señal emitida y señal reflejada, y a partir de ese desfase se obtiene el valor de la impedancia. Puede esperarse que los sensores tengan diferentes funciones de calibración en función de su frecuencia de operación; dicha frecuencia es el principal factor que afecta a las medidas de una sonda de capacitancia, provocando sensibilidad frente a cambios en las mencionadas condiciones del suelo. Según Robinson y otros [2008], estas sensibilidades se reducen al mínimo cuando la frecuencia es superior a los 500 MHz. Desde ese punto de vista, los sensores TDR, que actúan a frecuencias más altas (Walker y otros [2004], Chen y Or [2006]) han demostrado mayor precisión; los FDR, al operar en una frecuencia menor, son más sensibles a perturbaciones del entorno, como por ejemplo la temperatura, la textura o la conductividad eléctrica (Topp y otros [2000]; Zhang y otros [2004]; Bogena y otros [2007]). Por lo tanto, lo que los sensores TDR parecen a priori la mejor elección. 18 Capítulo 2. Calibración por Redes Neuronales Artificiales Sin embargo, aunque las sondas de capacitancia son mucho menos fiables, a cambio ofrecen una ventaja significativa: un coste inferior. Esto permite el uso de un número mayor de instrumentos, lo que implica redes mucho más densas (Mittelbach y otros [2011]). Dada la gran variabilidad espacio-temporal de la humedad en suelo, se puede postular que tener un gran número de sensores menos precisos pero más baratos pueden reducir el error de muestreo debido a variabilidad espacial (Teuling y otros [2006]). Los sensores capacitivos Decagon (Decagon Devices, Inc., Pullman, WA., EE.UU.) han tenido un relativo éxito en el mercado; sin embargo, sucesivos trabajos de evaluación y calibración en laboratorio (Blonquist y otros [2005], Jones y otros [2005], Bogena y otros [2007], Bandaranayake y otros [2007], Sakaki y otros [2008], Kizito y otros [2008]) han mostrado una gran sensibilidad al tipo de suelo, tensión de alimentación y conductividad eléctrica del suelo. La figura 2.3 muestra tres de estos sensores pertenecientes a la familia ECH2O: EC-20, 10HS y 5TE. Figura 2.3: Sensores de humedad de Decagon Devices, Inc. De arriba a abajo: EC20,10HS y 5TE. Modelado, simulación y control del uso del agua en la agricultura 19 2.1.2. Redes Neuronales Artificiales Las Redes Neuronales Artificiales (RNA) son un modelo matemático estocástico que se inspira en las redes neuronales naturales (biológicas), esto es, pretenden modelar el funcionamiento del sistema nervioso animal. Los primeros intentos de crear este tipo de herramienta comienzan con McCulloch y Pitts [1943], en los que se demuestra que un modelo simple con interruptores electrónicos es capaz de modelar con éxito una gran variedad de funciones analíticas. La figura 2.4 muestra el esquema básico de una red neuronal artificial. Figura 2.4: Esquema básico de una red neuronal con n entradas, una salida y una capa oculta de 30 perceptrones. Una RNA consta de un cierto número de nodos denominados «perceptrones» agrupados en dos o más capas comunicadas entre sí, de modo que cada nodo de una capa recibe la salida de todos los nodos de la capa anterior; en el caso de la figura, se parte de dos capas (la oculta, con treinta nodos, y la de salida) si bien el número de capas y de los nodos que cada una de ellas contiene puede ser ajustado de acuerdo con las necesidades del problema. Se parte de un primer nodo al que se suministra los datos entrada de un modelo; además, existe otro nodo a la salida que proporciona el resultado de la operación matemática que realiza la red; cada perceptrón lleva a cabo una función matemática cuya variable de entrada está formada por la suma ponderada (la ponderación es el peso W, correspondiente a cada uno de los nodos de la capa anterior) con un cierto umbral bde cada una de las salidas de los nodos de la capa anterior. Al comparar los datos de salida con el resultado esperado (la referencia), se produce una realimentación de los perceptrones, que ajustan sus pesos para aproximarse mejor al valor de referencia; a este proceso se le llama entrenamiento. Tras el entrenamiento, se llevan a cabo las fases de validación del modelo y, finalmente, de prueba. Al acabar, se dispone de una red completamente 20 Capítulo 2. Calibración por Redes Neuronales Artificiales caracterizada por sus pesos y umbrales que actúa como un modelo del proceso cuyas entradas y referencias se utilizaron para su creación. Con la evolución de los computadores electrónicos, matemáticos y neuropsicólogos avanzaron en la materia hasta que en 1958 Rosenblatt y Wightman (Rosenblatt [1958]) desarrollaron en el Instituto Tecnológico de Massachussets la primera neurocomputadora viable, el «Perceptron Mark 1». Minsky y Papert [1969] mostraron varias carencias en el modelo del perceptrón que detuvieron el desarrollo de las RNA hasta que Rumelhart y otros [1986] desarrollaron el método de propagación de errores hacia atrás (error back-propagation). Desde ese momento, como indican Zou y otros [2010], las RNA han experimentado un desarrollo expansivo en multitud de ámbitos gracias a su simplicidad, rapidez y rendimiento; a esta evolución ha contribuido decisivamente el continuo desarrollo de la electrónica y la informática. Hay que esperar hasta la segunda mitad de la década de 1990 para empezar a ver aplicaciones de RNA en los campos de agroecología e hidroecología, con el desarrollo de modelos de agroecosistemas de Schultz y Wieland [1997] o el sistema de regulación de fertilización de Broner y Comstock [1997]. Por su naturaleza, las RNA resultan particularmente interesante para el modelado de procesos no lineales, como por ejemplo la dinámica de la humedad y salinidad del suelo, estudiada por Morshed y Kaluarachchi [1998], o de la retención de agua y la conductividad eléctrica del suelo propuesta por Schaap y Leij [1998]. Más recientemente, Landeras y otros [2008] compararon el comportamiento de un modelo de evapotranspiración diaria basado en RNA con ecuaciones de modelado empíricas. Además, las RNA se han consolidado como una herramienta válida para modelar el sistema dinámico de la evolución de la humedad y salinidad del suelo, empleando para ello tanto RNA convencionales (Zou y otros [2010], Dai y otros [2011] como de alto orden Elshorbagy y Parasuraman [2008]. En este capítulo se persiguen dos objetivos: por un lado, proponer una mejora de la calibración en laboratorio mediante un enfoque mixto campo - laboratorio, exponiendo un protocolo de actuación a tal efecto. Por el otro, evaluar el rendimiento de redes neuronales artificiales en la calibración de sensores capacitivos en comparación con los métodos tradicionales para ajustes de curvas de calibración; la utilidad de las RNA para modelar procesos no lineales las convierte en un método más prometedor que el ajuste funcional y su simplicidad (relativa) de programación en un computador digital lo hace viable. Modelado, simulación y control del uso del agua en la agricultura 21 2.2. Material y métodos 2.2.1. Ámbito del ensayo Para llevar a cabo la calibración, se extrajeron cuatro columnas de suelo inalterado (105 mm de diámetro, 300 mm de altura) de puntos seleccionados de una parcela situada en la finca experimental «Tomejil», en el término municipal de Carmona (Sevilla). donde se comparan desde 1982 distintos sistemas de manejo de suelo. La figura 2.5 muestra una vista aérea de la parcela. Figura 2.5: Vista aérea de la parcela experimental en la finca «Tomejil». Dos de las columnas provenían de parcelas elementales bajo siembra directa y otras dos de parcelas bajo laboreo convencional. El suelo, formado sobre margas del mioceno, se clasifica como Haploxerept típico (U.S. DoA, Soil Survey Staff[1999]) con un contenido de arcilla de entre el 50% y el 64%. En dicha finca se ha instalado una red de sensores capacitivos Decagon; al ser la exactitud de esos sensores limitada, la necesidad de una calibración adecuada para éstos 28 Capítulo 2. Calibración por Redes Neuronales Artificiales RMSE =v t1 n n ∑ i=1 [yi−ˆ f(xi)] R2= n ∑ i=1[yi−ˆ f(xi)] n ∑ i=1(yi−¯ y) r2=R2 2−R2 1 1−R2 1 (2.2) Donde yison los valores de humedad obtenidos a partir de los datos gravimétricos, xilos valores raw de los sensores, ˆ f(xi) la salida del modelo y la media ¯ y=1 n n ∑ i=1yi. Además, R2 ison los coeficientes de determinación de dos modelos diferentes; r2, el índice de mejora de Nash y Sutcliffe, compara la bondad de dos ajustes diferentes informando de cuánto mejora el segundo con respecto del primero. Finalmente, se propone un modelo diferente: puesto que en cada columna hay tres sensores muy próximos en un espacio reducido, se puede asumir de forma razonable que sus medidas de humedad han de ser la misma. Desde ese punto de vista, puede ser interesante buscar una calibración múltiple: varios sensores para dar una sóla medida. Se diseñó una red neuronal artificial con tres entradas y una salida y se comparó con la referencia calculada a partir de los datos gravimétricos. 2.3. Resultados y discusión 2.3.1. Recopilación de datos La tabla 2.1 muestra los datos recogidos de las cuatro columnas de suelo. Modelado, simulación y control del uso del agua en la agricultura 29 Tabla 2.1: Datos de las columnas de suelo usadas; LC: laboreo convencional. SD: siembra directa. Columna LC-1 LC-2 SD-1 SD-2 Tara (g) 397 373 402 387 Peso saturado (g) 3780 3697 3539 3692 Vol. agua saturado (cm3) 1105 1101 992 1004 Volumen saturado (cm3) 2675 2596 2547 2688 Volumen seco (cm3) 2118 2196 1853 2034 Figura 2.11: Medidas de humedad y temperatura en una columna de Laboreo Convencional. Las figuras 2.11 y 2.12 muestran los resultados obtenidos de los diferentes sensores de humedad en suelo (medidas volumétricas) junto a las medidas gra- 30 Capítulo 2. Calibración por Redes Neuronales Artificiales vimétricas obtenidas y a las temperaturas tanto en la estufa como en el interior de dos de las columnas (una por manejo: la primera para Laboreo Convencional y la segunda para Siembra Directa). Puede comprobarse cómo la temperatura en el interior del cilindro sigue una dinámica de primer orden con respecto a la temperatura ambiente de la estufa y que conforme el cilindro se va secando, la constante de tiempo de esa dinámica se va reduciendo. Figura 2.12: Medidas de humedad y temperatura en columna de Siembra Directa. 2.3.2. Cálculo de las redes neuronales artificiales. Es relativamente fácil obtener una buena calibración para un sensor específico en unas circunstancias concretas. A título de ejemplo, la figura 2.13 muestra Modelado, simulación y control del uso del agua en la agricultura 31 cómo una calibración por redes neuronales puede dar lugar a un ajuste prácticamente perfecto en contraste con una calibración polinómica: en el caso de un sensor EC20, cuya calibración de fábrica sigue una norma lineal, una red de tan sólo diez neuronas mejora sustancialmente los resultados de una calibración polinómica cuadrática. Figura 2.13: Calibración específica para un único sensor; comparación con calibración polinómica. Las tablas 2.2 y 2.3 ilustran el efecto de tratar de emplear RNA (20 nodos) entrenadas con un sensor para otro sensor del mismo tipo en cada una de las cuatro columnas de suelo; se resaltaron las diagonales de las matrices de resultados, correspondientes a probar una RNA con el sensor específico para el que fue entrenada. Puede comprobarse una gran variabilidad en el comportamiento de los distintos modelos, observando que el coeficiente de determinación (con excepciones: véase el efecto de la red del 10HS del LC-1 con los datos del LC-2) resulta 32 Capítulo 2. Calibración por Redes Neuronales Artificiales Tabla 2.2: RMSE resultado de emplear una RNA entrenada con un sensor con los datos de los demás. Columna izquierda: sensor para el que se entrenó. Fila superior: columnas para las que se obtuvieron los datos. RMSE(m3/m3) Sensor Columna LC-1 LC-2 SD-1 SD-2 EC20 LC-1 0,00504 0,00504 0,00504 0,0560 0,0682 0,0574 LC-2 0,0442 0,00806 0,00806 0,00806 0,0281 0,0161 SD-1 0,0613 0,0412 0,00556 0,00556 0,00556 0,0207 SD-2 0,0436 0,0228 0,0182 0,00415 0,00415 0,00415 10HS LC-1 0,00503 0,00503 0,00503 0,0556 0,0682 0,0574 LC-2 0,0442 0,00806 0,00806 0,00806 0,0281 0,0161 SD-1 0,0613 0,0412 0,00560 0,00560 0,00560 0,0207 SD-2 0,0436 0,0228 0,0182 0,00415 0,00415 0,00415 5TE LC-1 0,00484 0,00484 0,00484 0,107 0,0775 0,0800 LC-2 0,0770 0,0143 0,0143 0,0143 0,0398 0,0417 SD-1 0,0676 0,0411 0,0112 0,0112 0,0112 0,0240 SD-2 0,0672 0,0428 0,0222 0,00965 0,00965 0,00965 Tabla 2.3: Coeficiente de determinación del modelo (análoga a 2.2). R2 Sensor Columna LC-1 LC-2 SD-1 SD-2 EC20 LC-1 0,998 0,998 0,998 0,682 0,428 0,490 LC-2 0,840 0,993 0,993 0,993 0,903 0,960 SD-1 0,693 0,828 0,996 0,996 0,996 0,933 SD-2 0,845 0,947 0,959 0,997 0,997 0,997 10HS LC-1 0,996 0,996 0,996 0,305 0,934 0,664 LC-2 0,611 0,977 0,977 0,977 0,530 0,375 SD-1 0,956 0,440 0,992 0,992 0,992 0,719 SD-2 0,856 0,314 0,912 0,995 0,995 0,995 5TE LC-1 0,998 0,998 0,998 0,917 0,955 0,968 LC-2 0,517 0,979 0,979 0,979 0,805 0,730 SD-1 0,627 0,829 0,985 0,985 0,985 0,911 SD-2 0,632 0,815 0,940 0,986 0,986 0,986 más alto cuando se utilizan datos de columnas de suelo del mismo manejo, especialmente para las de siembra directa. A la vista de estos resultados, se propuso calcular una única RNA para cada modelo de sensor: EC-20, 10HS y 5TE: uniendo las medidas de varios senso- Modelado, simulación y control del uso del agua en la agricultura 33 res se podría obtener una calibración que sirva en distintas condiciones para distintos sensores (y cuanto más diversas sean esas condiciones mejor será esa calibración general); a cambio, se pierde la precisión que proporciona la calibración específica. Con el propósito de calcular las RNA que se correspondan con los tres modelos de sensores, se configuró la red (en los tres casos) para que tuviera una única capa oculta de 30 nodos y se empleó un 70% de las muestras para el entrenamiento y un 15% para las fases de validación y pruebas. Para el cálculo de la red óptima se buscó la minimización del RMSE mediante el algoritmo de Levenbergh-Marquardt (Marquardt [1963]). La tabla 2.4 muestra los errores cuadráticos medios de las tres redes y los coeficientes de determinación correspondientes. En ella puede apreciarse un coeficiente de determinación que supera el 90% en los tres casos, llegando a alcanzar el valor de 0,96 en el caso del sensor EC20. Tabla 2.4: Resultados del entrenamiento de las RNA agrupando los datos por tipo de sensor. Sensor MSE R2 EC20 4,3679·10−40,955180 10HS 1,01492·10−30,9027 5TE 9,6792·10−40,9028 2.3.3. Evaluación de resultados La figura 2.14 muestra las curvas de regresión para el caso del sensor EC20. Pueden extraerse curvas muy similares para los tres casos: almacenando los valores de pesos y umbrales, se pueden definir las tres redes y puede procederse a probarlas con los sensores de cada caso. Para poder evaluar la bondad de estos resultados, se llevó a cabo un ajuste polinómico de orden tres para cada uno de ellos, obteniendo las expresiones mostradas en las ecuaciones 2.3, 2.4 y 2.5 (donde ysensor indica el sensor al que se refiere la función de ajuste) y se procedió a comparar los resultados. 34 Capítulo 2. Calibración por Redes Neuronales Artificiales Figura 2.14: Rectas de regresión para la red neuronal artificial para el sensor EC20. yEC20 = 0,1867(x−1157 287,4)3+0,02032(x−1157 287,4)2 +0,05832x−1157 287,4+0,173 (2.3) Modelado, simulación y control del uso del agua en la agricultura 35 y10HS = 0,006488(x−1197 168 )3+0,03271(x−1197 168 )2 +0,08133x−1197 168 +0,1662 (2.4) y5T E = 0,0007493 ·10−4(x−736 407,1)3 −0,01085(x−736 407,1)2 +0,1036x−736 407,1+0,2075 (2.5) La tabla 2.5 muestra los resultados de comparar la calibración con RNA para cada tipo de sensor con la calibración polinómica; puede observarse cómo las RNA mejoran el resultado de los ajustes polinómicos en prácticamente todos los casos, llegando a una mejora sensible en casos concretos (5TE en el LC1, 10HS en el SD1) y siendo sólo mínimamente inferior para los casos del sensor EC-20 en siembra directa. Tabla 2.5: Resultado de emplear la calibración global de los sensores sobre los datos de cada sensor. El error RMSE se mide en m3/m3. Columna Sensor RNA Polinómica RMSE R2RMSE R2 Laboreo convencional 1 EC-20 0,0389 0,877 0,0411 0,862 10HS 0,0294 0,930 0,0325 0,914 5TE 0,0493 0,802 0,0592 0,714 Laboreo convencional 2 EC-20 0,0152 0,977 0,0192 0,963 10HS 0,0144 0,951 0,0150 0,947 5TE 0,0336 0,886 0,0358 0,870 Siembra directa 1 EC-20 0,0352 0,848 0,0339 0,858 10HS 0,0661 0,768 0,0729 0,717 5TE 0,0260 0,917 0,0283 0,902 Siembra directa 2 EC-20 0,0098 0,985 0,0096 0,986 10HS 0,0276 0,882 0,0315 0,846 5TE 0,0158 0,961 0,0193 0,942 La figura 2.15 muestra la comparación gráfica entre ambas calibraciones en dos casos concretos: el EC20 en el cilindro 1 de laboreo convencional frente al 5TE en el cilindro 2 de siembra directa. Con respecto a la distribución de humedad en la columna de suelo, se puede ver en la figura 2.16 cómo la variación de ésta en los cilindros marcadores (en la gráfica indicados como C1, C2 y C3) con el diámetro de acuerdo con las 36 Capítulo 2. Calibración por Redes Neuronales Artificiales Figura 2.15: Comparación gráfica de calibraciones. condiciones descritas en la sección anterior es lo bastante reducida como para asumir que no debería afectar a los resultados. Los sensores están insertados en puntos donde se mide el valor medio, que viene representado en la figura por las líneas horizontales. Los niveles medios de humedad gravimétrica cumplen que θg1> θg2> θg3, lo que es coherente con las fechas de extracción. Modelado, simulación y control del uso del agua en la agricultura 37 Figura 2.16: Distribución de humedad gravimétrica en columnas de suelo. 2.3.4. Modelo de tres sensores. Finalmente se diseñó una única RNA en similares condiciones a las anteriores: una capa oculta de 30 perceptrones, pero con la diferencia de ser una red multivariable (tres entradas, una salida); por ello la capa de entrada, que en origen tenía un nodo, ahora tendrá tres, y los nodos de la capa oculta tendrán tres pesos cada uno. En estas condiciones se obtiene un coeficiente de determinación R2= 0,991 y un error cuadrático medio RMSE = 0,012(m3/m3). Las rectas de regresión pueden verse en la figura 2.17. Los resultados para cada uno de los cuatro casos pueden verse en la tabla 2.6. Puede observarse que en los cuatro casos se ha logrado un error de calibración mínimo y el índice de mejora de la eficiencia de Nash-Sutcliffe con respecto a los modelos obtenidos anteriormente supera el 57% para los sensores en columnas de siembra directa (86% en media) y el 91% para los de laboreo convencional (97% en media). 44 Capítulo 3. Modelo del contenido de humedad en suelo 3.2. Material y métodos 3.2.1. El modelo del cubo Un primer intento por modelar la dinámica del agua en el suelo en forma de ecuación de continuidad se debe a Thornthwaite [1948]. Este modelo se ha ido refinando a lo largo del tiempo; por ejemplo, Alley [1984] planteó los distintos componentes del balance de aguas. El modelo a partir del cual se desarrolla este capítulo fue propuesto por Rodríguez-Iturbe y otros [1999]: consiste en un modelo de tipo «leaky bucket», una ecuación de continuidad en la que se establece un balance de aguas, estipulando que el incremento del contenido de humedad en suelo varía linealmente con la suma de una serie de aportaciones y pérdidas. nZr dS dt =I(S,t)−L(S)−T(S)−E(S) (3.1) donde: nes la porosidad del suelo. Zres la profundidad de la zona de las raíces. S(t) es el contenido de humedad en suelo. I(S,t) es la infiltración de agua en el suelo; se obtiene sumando la lluvia y el riego (de haberlo) y restando la intercepción de la vegetación. L(S) son las pérdidas por percolación. T(S) son las pérdidas debidas a la transpiración de la planta. E(S) son las pérdidas por evaporación. Infiltración En el modelo original se trabaja en grandes extensiones de terreno (a nivel de cuenca) sin regadío, por lo que todas las infiltraciones se deben a precipita- Modelado, simulación y control del uso del agua en la agricultura 45 ciones. Estas precipitaciones se modelan como procesos estocásticos, de modo que se modela la llegada de tormentas como una variable aleatoria de Poisson y la intensidad de dichas tormentas (medidas en mm, esto es, l/m2) como una variable aleatoria exponencial. A su vez, a las precipitaciones hay que restar la intercepción debida a vegetación y el volumen no saturado de la zona radical: I(S(t− i),ti)=m´ ın[P(ti),nZr(1 −S(t− i))] (3.2) donde P(ti) son las precipitaciones y t− ies el instante inmediatamente anterior. Percolación En Laio y otros [2001] se describe el comportamiento de la percolación de la zona de las raíces como la conductividad hidráulica no saturada, descrita en la forma L(S) = Ksat eβ(S−Sf c)−1 eβ(1−Sf c)−1(3.3) donde Ksat es la conductividad hidráulica saturada y βes un parámetro del suelo. La capacidad de campo Sf c se define como el punto por debajo del cual la percolación es despreciable con respecto al resto de pérdidas. 46 Capítulo 3. Modelo del contenido de humedad en suelo Evapotranspiración Las «pérdidas» por transpiración no son tales: interesa que la planta transpire en la medida de sus posibilidades. Su expresión se obtiene de: T(S) =            0, S ≤Sw; S−Sw S∗−SwTmax Sw< S ≤S∗; Tmax, S > S∗. (3.4) donde Tmax representa la máxima tasa de transpiración que marca la demanda atmosférica en la unidad de tiempo. S∗es el umbral de saturación a partir de la cual la planta ya no puede transpirar más. Swes el umbral inferior (humedad residual) por debajo del cual la planta no tiene suficiente humedad para llevar a cabo el proceso de transpiración. Análogamente, la evaporación de la zona radical se puede obtener a partir de: E(S) =            0, S ≤Sh; S−Sh S∗−ShEmax Sh< S ≤S∗; Emax, S > S∗. (3.5) donde Emax representa la máxima evaporación en la unidad de tiempo y Shes el umbral por debajo del cual la tierra retiene la humedad, impidiendo su evaporación. La figura 3.1 muestra la suma de los tres efectos de pérdida: pueden observarse cinco regiones. 1. En la primera región, S < Sh, la humedad del suelo es tan pequeña que éste la retiene, impidiendo evaporación y transpiración. 2. En la segunda región, muy pequeña (Sh< S < Sw), existe evaporación pero la planta no puede transpirar. 3. En la tercera región (Sw< S < S∗), la evaporación y la transpiración actúan de forma lineal sobre el contenido de agua en la zona de las raíces. Modelado, simulación y control del uso del agua en la agricultura 47 4. En la cuarta región (S∗< S < Sf c), se alcanza una saturación evapotranspirativa; ni la atmósfera ni la planta son capaces de expulsar la humedad del suelo más allá de un límite teórico, la evapotranspiración máxima. 5. En la quinta región, S > Sf c, las pérdidas por percolación predominan sobre el resto de los efectos. Figura 3.1: Suma de los efectos de las pérdidas 3.2.2. Modificaciones El modelo descrito se enuncia con el objeto de aplicarse a escala de cuenca. El carácter estocástico de las precipitaciones en dicho modelo permite predecir 48 Capítulo 3. Modelo del contenido de humedad en suelo la forma en la que se va a comportar la variable del sistema (contenido de humedad en suelo) a lo largo del tiempo, típicamente un año (una campaña). Sin embargo, no hay nada que impida reducir la escala de trabajo desde la cuenca hasta la parcela, lo que permitiría aplicar el modelo a escala agronómica; si bien los experimentos se han conducido en cultivos de secano, no hay ningún motivo real por el que no se pueda aplicar a regadíos, con lo que el estudio a nivel de parcela permitiría tener un modelo del sistema suelo-planta-atmósfera (en lo sucesivo, SPA) que permita el diseño de controladores de riego que mantengan la humedad dentro de unos márgenes óptimos. Más aún, incluso dentro de una parcela, es posible dividir ésta en sectores, ya que estos no tienen por qué tener el mismo contenido de agua en suelo debido a su diferente orografía. Esta subsección detalla las modificaciones necesarias para llevar a cabo este cambio de enfoque. Ángulos de inclinación Si dividimos la parcela en sectores iguales (por simplicidad) de forma rectangular, se puede aplicar el modelo descrito a cada uno de éstos por separado sin necesidad de alterarlo. Sin embargo, esto no será correcto debido a un factor adicional no tenido en cuenta originalmente: la orografía del terreno. Esto es, si dos sectores contiguos están a diferentes alturas, es de esperar que el agua obedezca a la ley de la gravedad y se transmita del sector más alto al más bajo, al menos en parte. Es por ello que necesitamos añadir una variable adicional al sistema: la altura media del sector, h. Evidentemente, para una precisión óptima desde el punto de vista teórico tendríamos que trabajar con sectores de tamaño infinitesimal y alturas instantáneas, pero eso no resultaría práctico y ha de alcanzarse un compromiso respecto al tamaño razonable del sector. A menudo, lo más recomendable es aprovechar la geometría del sistema de regadío y tener tantos sectores en la parcela como sectores de riego diferentes el agricultor desee instalar. Una vez obtenida la altura del sector h, resulta fácil obtener el ángulo de inclinación entre dos sectores, θ: θ= arctan l h1−h2(3.6) Modelado, simulación y control del uso del agua en la agricultura 49 Donde les el lado del sector (distancia entre los puntos medios de ambos sectores) y hila altura de cada sector. Para cada sector se calculará el ángulo de inclinación con respecto a los cuatro sectores contiguos. Una alternativa a esta geometría, siguiendo a Toth [1972], sería emplear celdas hexagonales, la forma más eficiente de teselar un plano en el espacio euclídeo. Sin embargo, llevar a cabo ese tipo de división de parcela en campo añade excesiva complicidad al sistema. Flujo lateral Si tomamos un sector en tres dimensiones como un paralelepípedo de lados (a,b) y profundidad Zr, es posible definir un flujo de pérdidas separado en dos componentes:  FL=FLuˆ u+FLzˆ z(3.7) Donde ˆ ues un vector unitario en la dirección de inclinación (esto es, tangente a la superficie) y ˆ zes el vector unitario vertical. Asumiendo flujo constante, se podrán separar las pérdidas en tres componentes: aporte lateral en dirección x, aporte lateral en dirección y y pérdidas por percolación. Naturalmente, puesto que cada sector tiene cuatro sectores adyacentes, se dividirá a su vez el sector a su vez en cuatro componentes, cada uno de los cuales tendrá un aporte a o desde sectores adyacentes. Aporte lateral a los sectores en la dirección del eje x: Lx=L 4zr asen(θx) si θx>0 Aporte lateral a los sectores en la dirección del eje y: Ly=L 4zr bsen(θy) si θy>0 Pérdidas por percolación: Lz=L−Lx1−Lx2−Ly1−Ly2 Naturalmente, si el ángulo de inclinación es cero o negativo, los términos de aportes laterales se anulan. El siguiente paso es añadir a la infiltración los aportes de estos sectores contiguos, de modo que partiendo de la ecuación 3.2 50 Capítulo 3. Modelo del contenido de humedad en suelo la infiltración y las pérdidas quedan como se muestra en la ecuación 3.8, donde los idel sumatorio hacen referencia a los elementos contiguos y lies la longitud del lado que actúa como frontera con el elemento i. I(S(t−),t) = min{P(t),nZr[1 −S(t−)]}+ 1..,4 ∑ θi<0 Li Li=       L 4zr lisen(θi) si θi>0 0 siθi≤0 (3.8) Escorrentía Comoquiera que el flujo lateral es un componente de las pérdidas por percolación, por lo visto en la figura 3.1 su efecto será pequeño excepto en momentos de saturación, que en principio un control adecuado trataría de evitar (o sobre los que, en el peor de los casos, no actuaría): «el flujo lateral a través del suelo es despreciable en entornos semi-áridos», según Bull y Kirkby [2002]. Sin embargo, hay otro efecto de la diferencia de alturas que sí resulta relevante: la escorrentía (Reaney y otros [2007]), esto es, el flujo lateral de agua sobre la superficie del terreno, sin que se produzca infiltración en éste. La escorrentía, como es lógico, sólo se producirá cuando exista precipitación o bien cuando el terreno esté tan saturado de agua que no sea capaz de filtrar el aporte en el ciclo de tiempo anterior. Dependerá, naturalmente, del tiempo y del ángulo de inclinación. Para simplificar el problema, se asumirá que la dependencia temporal se halla implícita en las precipitaciones P(t) y en la evolución del contenido de agua en suelo S(t). Asumiendo una dinámica de primer orden, caracterizamos la escorrentía como: Esci(t) = I(1 −e−S(t)/τesc)θi ∑4 j=1,θj>0θj (3.9) donde τesc es el coeficiente de escorrentía, una constante que habrá de añadirse al modelo y obtener mediante el proceso de calibración. Modelado, simulación y control del uso del agua en la agricultura 51 Lo que nos lleva a un problema circular: la escorrentía depende de la infiltración, pero a su vez la infiltración depende de la escorrentía. Para resolver este conflicto, la escorrentía se calculará a partir del valor de la infiltración en el instante anterior t−, quedando la expresión de la ecuación 3.10 donde los sumatorios se refieren a los sectores contiguos. I(S(t−),t)=m´ ın{P(t),nZr[1 −S(t−)]}+ 1..,4 ∑ θi<0 Li+ 1..,4 ∑ θi<0 Esci(t)− 1..,4 ∑ θi>0 Esci(t) (3.10) 3.2.3. Entorno de simulación Los experimentos se llevaron a cabo a lo largo de dos años en una parcela perteneciente a en la finca «Tomejil», situada en el término municipal de Carmona (Sevilla). La parcela se divide en dieciocho sectores tal y como muestra el diagrama de la figura 3.2. Figura 3.2: Separación en sectores de la parcela de estudio El modelo descrito en 3.2.2 se programó en entorno Matlab (©The Mathworks, inc. Natick, MA, EE.UU.). El programa consta de cuatro partes principales: 52 Capítulo 3. Modelo del contenido de humedad en suelo 1. La célula. Es un objeto que modela un sector de la parcela y sus propiedades son sus dimensiones (xey), su altura sobre el nivel del mar (las alturas en metros se pueden observar en la figura 3.2) y los ángulos de inclinación θcon respecto a sus sectores vecinos. 2. Las tablas. La interfaz del programa con el usuario, donde puede verse el estado de la simulación en todo momento. Existe una tabla donde se describen las posiciones de las células en el plano (la parcela), otra donde encontramos los aportes laterales y de escorrentía para cada sector en el instante de tiempo de la simulación y -finalmenteuna tabla adicional que recoge el contenido de humedad en suelo (salida del modelo y variable de interés) en el instante de tiempo de simulación. 3. La función f. El modelo descrito se puede resumir en dS dt =f(S,t), así que es imprescindible una función que calcule cuánto vale esa f(S,t) en cada caso. Devolverá las dos tablas de información para un instante de tiempo: el contenido de humedad en suelo y los términos de escorrentía y aporte lateral que se intercambian entre los distintos sectores. 4. El algoritmo de simulación. Tomando como entrada las precipitaciones en el intervalo de tiempo T, el algoritmo lleva a cabo un procedimiento de integración numérica empleando la herramienta ode45 de Matlab®; dicha herramienta utiliza una fórmula explícita de Runge-Kutta 4,5 (de ahí su nombre: Ordinary Differential Equation 4,5) conocida como la pareja Dormand-Prince (Dormand y Prince [1980]). Como en este sistema las lluvias actúan como entrada y no son predecibles, la función es llamada una vez por cada intervalo de integración (tiempo de muestreo) Tm. De este modo, no sólo se obtiene la humedad en cualquier instante, sino que se puede obtener también otra información (por ejemplo, la escorrentía) para cada paso de integración. El programa completo crea las tablas, ejecuta el algoritmo y representa a lo largo del tiempo el contenido de humedad en suelo para cada sector y la media de la parcela. 3.2.4. Calibración Para calcular los parámetros del suelo, se aprovechó un experimento previo, desarrollado a lo largo de los años 2003 y 2004 por Jiménez-Bocanegra y Modelado, simulación y control del uso del agua en la agricultura 53 otros [2005] en una parcela adjunta, de idénticas características de suelo y cultivo aunque diferentes dimensiones. El motivo del uso de dicho experimento es que en dichas parcelas se hayan instaladas sondas de tipo EnviroScan (SENTEK Technologies Ltd., Stepney, Australia) de mayor precisión y fiabilidad (Arauzo y otros [2007]). Esta parcela adjunta consta de cuatro sectores de dimensiones 75 x 15 m, situados a las siguientes altitudes sobre el nivel del mar tal y como se muestra en la tabla 3.1. Tabla 3.1: Esquema de la parcela adjunta. Las altitudes se miden en metros. 80 80.5 81 81.5 Los datos de precipitación se obtienen a partir de una estación meteorológica perteneciente a la Red de Información Agroclimática del IFAPA y puesta a disposición del público. Dicha estación se halla ubicada a menos de cien metros de la parcela. Esos datos se utilizaron como entrada del sistema; los datos de humedad relativa obtenidas por las sondas EnviroScan se comparan con la salida del sistema y se ajustan los parámetros mediante el algoritmo de punto interior del Toolbox de optimización de Matlab®(Mehrotra [1992]). 3.3. Resultados Una vez se han aplicado todos los cambios indicados en la ec. 3.10 al modelo original de la ec. 3.1, se obtiene la expresión de la ec. 3.11. Hecho esto, se puede proceder a identificar los parámetros del modelo. nZr dS dt = m´ ın{P(t),nZr[1 −S(t−)]}+ 1..,4 ∑ θi<0 Li+ + 1..,4 ∑ θi<0 Esci(t)− 1..,4 ∑ θi>0 Esci(t)−L(S)−T(S)−E(S) (3.11) 60 Capítulo 4. Control del contenido de agua en suelo Bleby [2006]). Este proceso de redistribución del agua en el suelo por las raíces de las plantas adquiere una importancia relevante en el campo de la hidroecología, a medida que se va avanzando en su comprensión hasta el punto de que se le reconoce una clara influencia en el riego de los cultivos asociado a situaciones de clara heterogeneidad de los perfiles de humedad del suelo (Ryel y otros [2002]). Una eficiente observación de la humedad del suelo no tiene interés únicamente para el estudio de la interacción entre el clima, el suelo y la vegetación, sino que también supone un avance significativo en aplicaciones agroambientales que requieren datos sobre su dinámica espacio-temporal. Un ejemplo de ello es la comparación de la dinámica de la humedad del suelo en experimentos con distintos sistemas de manejo (Vanderlinden y otros [2008]) donde resulta técnicamente inviable obtener resultados fiables y concluyentes debido a varios factores: la variabilidad espacial de las propiedades del suelo. la gran incertidumbre asociada a las observaciones realizadas con sensores electromagnéticos (Hignett y Evett, 2008). limitaciones de carácter tecnológico y económico en la implementación de estos sistemas de medida (limitación de la longitud de los cables, impedimentos de las labores del suelo en los componentes superficiales, el elevado coste del almacenamiento de los datos in situ debido al gran número de unidades requeridas, etc.). la imposibilidad de realizar periódicamente observaciones destructivas simultáneamente en un gran número de puntos como alternativa al uso de sensores de humedad. Otro ejemplo es la creciente tendencia de apoyar decisiones sobre la programación del riego en datos de humedad del suelo (Thompson y otros [2007]), tanto en cultivos hortícolas como en cultivos herbáceos o leñosos, con el fin de optimizar al máximo el uso del agua de riego y fertilizantes en situaciones, cada vez más frecuentes, de suministro limitado de agua (Fereres y Soriano [2007]). La capacidad de registrar continuamente la distribución espacial de humedad del suelo a escala de parcela en tiempo real también abre la puerta al riego de precisión (Sadler y otros [2005], DeJonge y otros [2007]), aplicando agua de manera diferenciada en zonas de manejo o mediante equipos que permiten la aplicación de un riego localizado de alta precisión a escala de parcela (Al-Kufaishi y otros [2006]). Modelado, simulación y control del uso del agua en la agricultura 61 En este sentido, existen multitud de modelos a diferentes escalas para el modelado del agua en suelo. En Xu [1999] se propone tomar modelos a gran escala y hacer un re-escalado para aplicarlos a escalas más pequeñas; en este trabajo se toma este principio para, a partir de un modelo estocástico tradicionalmente empleado a escala de cuenca, modelar el contenido de agua en suelo a escala de parcela o incluso de sectores en los que se divida una parcela. Se consigue así un modelo de parámetros distribuidos, tal y como se plantea en Guo y otros [2002], que permite estudiar la interacción entre distintas secciones de una parcela. Como respuesta a las dificultades técnicas mencionadas aparece el uso de las redes inalámbricas de sensores que permiten el estudio de la humedad en un número de puntos aislados entre sí como representación del estado del conjunto de una parcela. El empleo de dicha tecnología no es una idea nueva en agronomía (Akyildiz y otros [2002]), pero en muchos aspectos se puede considerar que se encuentra en su fase inicial de utilización (Srivastava [2010]). En la actualidad, la mayoría de las redes inalámbricas de sensores instaladas pueden agruparse en dos grupos según el estándar de comunicación que empleen: WLAN (Wireless Local Area Network, IEEE 802.11) o WPAN (Wireless Personal Area Network, IEEE 802.15), que a su vez puede separarse a grandes rasgos en redes basadas en tecnologías Bluetooth y ZigBee. Esta última es la tecnología que parece más prometedora debido a su bajo consumo de energía y la sencillez de su configuración de red. Su principal inconveniente es que la normalización en el mercado no está completa, lo que da lugar a multitud de variantes propietarias del estándar (Wang y otros [2006]). A pesar de que la instalación de redes inalámbricas presenta dificultades y riesgos, su aplicación a las explotaciones agrícolas puede servir como punta de lanza para el empleo del uso de nuevas tecnologías en la agricultura. Una vez que ya se dispone de una infraestructura funcional para la recogida de información y un modelo matemático para procesar esa información, el siguiente paso es llevar a cabo una planificación del riego que optimice la productividad de agua (Arriaga y otros [2012]). La teoría de control (Ogata [2001]) se ocupa del diseño de controladores automáticos y reguladores que permitan obtener un mejor rendimiento de los sistemas dinámicos, incrementando la productividad y reduciendo tanto la carga de trabajo como los costes asociados. Arbat y otros [2013] presenta un software para simular patrones de penetración del agua en el suelo a partir de riego por goteo. Precisamente, ésta es una de las carencias de la investigación agronómica actual, ya que, a pesar de que se han realizado numerosos esfuerzos para modelar el continuo SPA, hay escasa literatura sobre control en estos sistemas. La teoría 62 Capítulo 4. Control del contenido de agua en suelo de control, tan extensa y exitosamente aplicada en otras áreas de la ciencia y la industria, apenas ha sido introducida en ámbitos agronómicos. El objetivo de este capítulo es, partiendo del modelo obtenido en el capítulo anterior, aplicar técnicas clásicas de control automático para lograr una regulación automática del contenido de agua en suelo, calculando controladores mediante aproximación lineal y probándolos en el entorno de simulación adecuado. 4.2. Material y métodos 4.2.1. Determinación del punto de operación El modelo no lineal del sistema, tal y como se presentó en el capítulo 3, se muestra en la ecuación 4.1. nZr dS dt = m´ ın{P(t),nZr[1 −S(t−)]}+ 1..,4 ∑ θi<0 Li+ + 1..,4 ∑ θi<0 Esci(t)− 1..,4 ∑ θi>0 Esci(t)−L(S)−T(S)−E(S) (4.1) A la hora de llevar a cabo el diseño de un controlador automático para un sistema no lineal como éste, existen varias técnicas Ogata [1996]. Sin embargo, la más sencilla consiste en linealizar el sistema: esto es, crear un modelo lineal que aproxime el sistema en torno al punto de operación; de este modo, es posible diseñar controladores pensados para sistemas lineales, más simples de diseñar y sintonizar. Ésta es la aproximación que se ha escogido a lo largo de este capítulo. El primer paso para llevar a cabo esa linealización es determinar el punto de operación o funcionamiento. Para una primera aproximación, se ignoraron los efectos de escorrentía y laterales y se trabajó tan sólo con el modelo del cubo inicial 3.1. El objetivo es Modelado, simulación y control del uso del agua en la agricultura 63 que la humedad del suelo S sea superior a SW, lo que implica que la planta pueda transpirar, y al mismo tiempo inferior a Sf c para que no se produzcan pérdidas (y, por lo tanto, tampoco se produzcan aportes de flujo lateral). A la vista de los valores obtenidos en la calibración, vistos en la figura 3.4, se buscará un contenido de humedad en suelo de 0,45. Sólo al concluir el diseño del controlador, se simula su rendimiento con el modelo completo de parámetros distribuidos, según plantea la ecuación 4.1. En el modelo se distinguen parámetros para las dos estaciones, la húmeda y la seca. El procedimiento a seguir es idéntico para ambos conjuntos de parámetros, por lo que sólo se detallan los cálculos de la estación húmeda. También por simplicidad, la profundidad de la zona radical, que se distingue para los dos tipos de manejo de suelo (Laboreo Convencional / Siembra directa), se toma como la media de ambos valores, Zr= 82,5cm. A partir de este punto de polarización, se procede a la linealización del sistema, obteniendo una función de transferencia: Y(s) U(s)=G(s) (4.2) Donde U(s) es la entrada del sistema, que se corresponderá con la infiltración I, e Y(s) la salida, que se corresponderá con el contenido de humedad en suelo S. 4.2.2. Diseño de controladores Aunque se ha definido el sistema como continuo, necesariamente habrá que trabajar con un sistema discreto puesto que no se tiene acceso a datos de precipitación continuos. Por ello, el sistema controlador será diseñado como discreto con tiempo de muestreo T= 1u.t. como se muestra en la figura 4.1. Puede observarse en dicha figura cómo se deshace el cambio de variable para que la salida se corresponda con la variable física que deseamos controlar, el contenido de humedad en suelo. 64 Capítulo 4. Control del contenido de agua en suelo Figura 4.1: Diagrama de control discreto del sistema linealizado Control Proporcional El control proporcional (P) es el caso más simple de control por realimentación de errores en la salida. Dada la función de transferencia G(s) del sistema linealizado (ecuación 4.2), se puede definir un C(s) = Pque estabilice la salida Y(s) dada una referencia R(s), como muestra la ecuación 4.3. Y(s) = G(s)C(s)(R(s)−Y(s)) Y(s)[1+ C(s)G(s)] = C(s)G(s)R(s) Y(s) = CG(s) 1+ CG(s)R(s) (4.3) La transformada de Laplace de una entrada en escalón de amplitud r, r(t) = r,t > 0, es R(s) = r/s. Aplicando el teorema del valor final, que establece que l´ ım t→∞y(t)=l´ ım s→0sY (s) se obtiene la expresión de la ecuación 4.4. Modelado, simulación y control del uso del agua en la agricultura 65 yr.p. = l´ ım s→0rP /(as) 1+ P /(as)= l´ ım s→0 rP P+as =r(4.4) La salida en régimen permanente coincide con la referencia. Esto era previsible, dado que el sistema linealizado es un integrador, lo que debe eliminar el error en régimen permanente frente a un escalón. Control Proporcional Integral (PI) El controlador PI añade, al término proporcional, un término integral que actúa eliminando -en el caso de que existael error en régimen permanente frente al escalón. Al ser el sistema un integrador, este efecto se acentúa eliminando el error en régimen permanente frente a una rampa. La ecuación en el dominio de la frecuencia de un controlador PI se muestra en la ecuación dada por la ecuación 4.5. C(s) = Kp Tis+1 Tis(4.5) Como el sistema es discreto, muestreado cada T= 1/96 días, se aplica una discretización de Euler hacia adelante (Åström y Hägglund [2006]), resultando la expresión de la ecuación 4.6. s=z−1 T C(z) = U(z)/E(z) = Kp Tiz+T−Ti Tiz−Ti =Kp Ti+(T−Ti)z−1 Ti−Tiz−1 (4.6) 66 Capítulo 4. Control del contenido de agua en suelo Llevándolo al dominio del tiempo mediante la aplicación de la transformada inversa Z, se obtiene la ecuación en diferencias 4.7. u(k) = u(k−1)+ Kpe(k)+(T /Ti−1)Kpe(k−1) (4.7) Control Proporcional Integral Derivativo (PID) Para conseguir que una respuesta más rápida se introduce un término derivativo, formando un controlador PID. Análogamente al caso anterior, la ecuación de un PID discreto por Euler hacia adelante en forma paralela viene dada por Åström y Hägglund [2006], que se muestra en la ecuación 4.8. C(z) = U(z)/E(z) = Kp[1+ T Ti(z−1) +Td(z−1) T+αTd(z−1)](4.8) Operando, se obtiene la expresión de la ec. 4.9 (1−z−1)U(z) E(z)=Kp(1−z−1)+KpT Ti z−1+KpTdP Tz−1 1−z−1+αTd ... T z−1(1−z−1)U(z)+ αTd(1−z−1)2U(z) =KpE(z)[Td(1 −z−1)+ +Tdz−1(1−z−1) +αTd(1 −z−1)2+ +T2 Ti z−2+αTd/Tiz−1(1−z−1)] (4.9) Para pasar al dominio del tiempo y conseguir la expresión en forma de ecuación en diferencias, será necesario hacer la transformada inversa Z, para lo que se aplican ciertas propiedades de dicha antitransformada mostradas en la ecuación 4.10. Modelado, simulación y control del uso del agua en la agricultura 67 X(z)(1−z−1)Z−1 −−−→x(t)−x(t−1) X(z)(1−z−1)2Z−1 −−−→x(t)+ x(t−2) X(z)z−1Z−1 −−−→x(t−1) X(z)z−2Z−1 −−−→x(t−2) (4.10) Aplicando estas fórmulas de transformada inversa, la expresión en el dominio del tiempo queda: u(k) = 1 αTd [−T u(k−1)+ (T−αTd)u(k−2)+ KpTd(1+α)e(k)+ +Kp(αTd Ti−Td+T)e(k−1)+ (KpTdα+Kp T2 Ti−KpαTd Ti−KpT)e(k−2)] (4.11) 4.3. Resultados 4.3.1. Linealización del sistema En el punto de operación escogido S0, como se puede ver en la figura 3.4 y muestra la ecuación 4.12., la evapotranspiración y las pérdidas por percolación permanecen constantes. ET (S0) = E(S0)+T(S0) = Emax +Tmax = 0,3 L(S0)≈Ksat eβ(S0−Sf c)−1 eβ(1−Sf c)−1= 0 (4.12) Asimismo, I(S,t) = min[P(t),nZr(1−S(t−))]. En el experimento de calibración realizado, S(t)<0,75, ∀t, por lo que se usará como cota máxima (y, por lo tanto, 68 Capítulo 4. Control del contenido de agua en suelo 0,25 como cota mínima). Dados los valores obtenidos para nZr en la calibración, I(S,t) = P(t) si P(t)<5,44cm. Por lo tanto, es posible eliminar la dependencia de S(t) en I(S,t)→I(t). Llevando estos valores a la ecuación 3.1, se obtiene: dS dt = [I(t)−ET (S)−L(S)]/nZr=I(t)−0,3 nZr Si aplicamos la transformada de Laplace para pasar al dominio frecuencial, el término −0,3 nZr presenta un problema: genera una dinámica no controlable. Por ello, se llevará a cabo un cambio de variable de entrada definiendo la variable auxiliar x(t) = I(t)−0,3 que será la nueva señal de entrada. Para evitar confusiones entre la variable objetivo S(t), contenido de humedad en suelo, y la variable de la transformación de laplace, s, se aplica el cambio de nomenclatura y(t) = S(t). Estos cambios se muestran en la ecuación 4.13. y(t) = S(t) x(t) = I(t)−0,3(4.13) Con lo que la ecuación diferencial del modelo resulta: dy(t) dt =x(t) nZr =x(t) 82,5(4.14) Pasando al espacio de Laplace, se obtiene: Y(s) X(s)=G(s) = 1 82,5s(4.15) Para evaluar la bondad de nuestro modelo, comparamos el comportamiento del sistema linealizado y el sistema completo, para un sólo sector, con varias referencias de entrada. La primera prueba es para la zona lineal: comprobar cómo entre contenidos de agua en el suelo de entre el 40% y el 55%, esto es, evapotranspiración máxima sin influencia de las pérdidas por percolación, am- Modelado, simulación y control del uso del agua en la agricultura 69 bos modelos coinciden. La figura 4.2 muestra el sistema en dos circunstancias opuestas: a la izquierda, cómo partiendo de una condición inicial de 0,4m3/m3, la salida crece al serle aplicado una entrada constante de 5mm/dia; a la derecha, el proceso opuesto al partir de una condición inicial de 0,55m3/m3sin precipitación. La linealización en esta región responde, como era de esperar, al comportamiento del modelo completo. Figura 4.2: Comportamiento de modelos (no lineal y linealizado) frente a una señal constante de 5mm/dia (izquierda) y 0 (derecha). Para comprobar cómo el modelo funciona fuera de esta zona de seguridad, se programó una segunda simulación, en este caso planteada a lo largo de un periodo de 150 días; partiendo de una condición inicial de 0,5m3/m3, se aplica un escalón unitario constante en la entrada entre los días 50 y 85. Puede observarse en la figura 4.3 cómo el comportamiento fuera de la zona lineal, si bien no se corresponde exactamente con el del sistema completo, sí que supone una aproximación satisfactoria. 76 Capítulo 4. Control del contenido de agua en suelo se alcance la ganancia última Ku, en la que la salida oscila con amplitud constante; sistema estable pero no asintóticamente estable. Una vez obtenida dicha ganancia y el periodo de oscilación de la salida para ella, Tu, se fijan las ganancias para los términos proporcional, derivativo e integral. El control ZieglerNichols clásico es un control de gran robustez, pero en este caso y a la vista de la inestabilidad del primer controlador se busca un control aún más robusto. De acuerdo con McCormack y Godfrey [1998], se escogen unas ganancias que eliminen la sobreoscilación: Kp= 0,2Ku;Ki= 2K∗p/Tu;Kd=KpTu/3. Figura 4.10: Simulación del sistema con un control PID sobre los tres primeros sectores La arquitectura de este PID será también más simple, a partir de Ogata [1996] la expresión de un PID discreto se muestra en la ecuación 4.16. u(t) = u(t−1)+Kp[(Kp+KiT+Kd T)e(t)+ (−Kp−2Kd T)e(t−1)+ Kd Te(t−2)](4.16) Modelado, simulación y control del uso del agua en la agricultura 77 En esas condiciones, se simula el sistema y se obtienen los resultados mostrados en la figura 4.10. Puede observarse cómo este controlador PID es más rápido que el PI, elimina el error en régimen permanente y prácticamente elimina la sobreoscilación, si bien su reacción frente a perturbaciones es un poco más elevada que en el control PI. 4.4. Conclusiones A partir del modelo calculado en el capítulo anterior, y del entorno de simulación construido al efecto, se han utilizado técnicas clásicas de control lineal para regular el contenido de agua en suelo, planteándose controladores individuales para cada celda del sistema. Se obtuvieron las dimensiones y altitudes en una parcela de dieciocho celdas situada en una finca experimental propiedad del IFAPA, y se simuló utilizando como entradas tanto la salida de los controladores como las precipitaciones obtenida a partir de una estación meteorológica cercana. El sistema se linealizó en torno a un punto considerado óptimo, resultando un sistema de primer orden. El sistema es estable, como cabía esperar: en ausencia de entradas, la humedad de suelo tiende a un valor próximo a cero (humedad residual). Se diseñaron tres controladores mediante ajuste computacional: los controladores proporcional (P), proporcional integral (PI) y proporcional integral derivativo (PID). En los dos primeros casos, se consiguió la estabilización del sistema a pesar de la presencia de perturbaciones -las precipitacionesque pueden llegar a ser muy fuertes; hay que tener en cuenta que por la propia naturaleza de la acción de control, ésta se halla limitada (saturada) a valores positivos: no se puede regar menos de cero. En el caso del controlador PID, la linealización del sistema no fue lo bastante buena como para resistir una acción derivativa, que inestabilizó el sistema; sin embargo, un ajuste empírico siguiendo la técnica clásica de Ziegler Nichols permitió resolver ese inconveniente estabilizando y regulando el sistema satisfactoriamente. La distinta altitud de las celdas de la parcela muestra cómo para el control proporcional, a pesar de la estabilización, el valor final de la humedad para un 78 Capítulo 4. Control del contenido de agua en suelo control proporcional varía según la celda (especialmente cuando ocurren precipitaciones). Por ello, el siguiente controlador a probar debería ser un controlador distribuido que para cada celda tenga en cuenta las aportaciones laterales de las celdas de alrededor. Sin embargo, los controladores PI y PID eliminan ese problema, aunque la reacción de éstos no es uniforme para todos los sectores, por lo que un controlador distribuido podría mejorar sus respuestas. El siguiente paso debería ser, evidentemente, probar estos controladores en un ensayo de campo real, instalando riego automático controlado en las celdas de la parcela donde se ha planteado. Parte II Modelado de los efectos del Riego Deficitario sobre cultivos leñosos 79 Modelado, simulación y control del uso del agua en la agricultura 81 Publicaciones Parte del contenido de los cuatro capítulos que constituyen esta segunda parte ha sido publicado en forma de artículos. García Tejero, I.F, Arriaga, J., Durán Zuazo, V.H., J. y Muriel Fernández, J.L. Enero 2013. «Predictiong crop-water production functions for long-term effects of deficit irrigation on citrus productivity (SW Spain).» Archives of Agronomy and Soil Science, Vol. 59, Is. 12, pp. 1591-1606. DOI: 10.1080/03650340.2012.748983 García Tejero, I.F., Durán Zuazo, V.H., Arriaga, J., Muriel Fernández, J.L. Julio 2012. «Impact of water stress on citrus yield.» Agronomy for Sustainable Development, Vol. 32, Is. 03, pp. 651-548. DOI: 10.1007/s13593-011-0060-y. Arriaga, J., García-Tejero, I.F., Muriel-Fernández, J.L., Durán-Zuazo, V.H. and Rubio, F.R. 2014. «Modeling, simulation and control of irrigation on young almond trees.» Acta Horticulturae (ISHS), vol. 1038, pp. 479-486. 82 Introducción La agricultura representa aproximadamente el 70% del total de agua dulce consumida en todo el mundo (García-Tejero y otros [2011a]), cifra que asciende a un 75% en el caso de España según el Ministerio de Medio Ambiente, Rural y Marino (Horizonte 2015 [2010]). La superficie de tierra dedicada al cultivo de regadío se ha incrementado en nuestro país desde los 3,2 millones de hectáreas de 1994 hasta los 3,5 de 2014 (Ministerio de Agricultura, Alimentación y Medio Ambiente ESYRCE [2014]), de los cuales un 29,26% se hallan en Andalucía. El agua destinado a regadío ascendía en 2008 según el citado Horizonte 2015 [2010], a 15.313 hm3al año, un 10% de ellos dedicado al cultivo del cítrico (MARM [2008]). De acuerdo con Margat y Vallée [1997], los recursos de agua se están volviendo cada vez más escasos; este fenómeno parece particularmente acusado en regiones del sur del Mediterráneo. Históricamente los agricultores han venido empleando grandes volúmenes de agua, excediendo ampliamente las necesidades de los cultivos tanto perennes como anuales (Shideed y otros [2005]), produciéndose sobrerriegos de entre el 30% y el 50% (Hamdy y Katerki [2006]). Mientras tanto, el contexto de cambio climático global amenaza con provocar un serio impacto en los recursos hídricos y en la agricultura en los próximos años. Es por eso que varios autores han tratado de evaluar las posibles consecuencias de este cambio en diversas regiones del mundo en los últimos años (Olesen y Bindi [2002], Aggarwal [2003], Jones y Thornton [2003]). Los modelos de cambio climático son múltiples y a veces incluso contradictorios, pero existe un cierto consenso en ellos en cuanto a predecir un aumento de la temperatura media global en el próximo siglo, incremento que según el modelo oscila entre los 1,1◦Cy los 6,4◦Centre los años 1990 y 2100 (Intergovernmental Panel for the Climate Change, Solomon y otros [2007]) y que vendrá acompañado de cambios en los patrones espacio-temporales de precipitaciones (Southworth y otros [2000], Räisänen [2001]). Los efectos de ese cambio serán muy probablemente más acusados en regiones áridas y semiáridas, que ya se encuentran en un delicado equilibrio que los hace más sensibles a cualquier cambio climático, y en las que los recursos hídricos ya son escasos. Incluso regiones semi-áridas que gracias a sus mayores desarrollo socioeconómico y acceso a tecnología, o incluso por disponer de climas relativamente más lluviosos, parecen disponer de suministros suficientes de recursos hídricos para destinar a la agricultura, es muy probable que sufran Modelado, simulación y control del uso del agua en la agricultura 83 los efectos adversos del mayor consumo de agua asociado a un clima más cálido IPCC2001b [2001]. Es por ello que el uso sostenible de los recursos hídricos en la Europa mediterránea es crucial para mitigar muchos de estos problemas de escasez y mal uso del agua, especialmente en áreas del sur de España (Zacharias y Koussouris [2000], Mariolakos [2007], Downward y Taylor [2007]). Al no verse limitados por la falta de disponibilidad de luz (antes bien, se dispone de sobreabundancia), la disponibilidad de agua es el principal factor limitante para el desarrollo de los cultivos. Es por eso que resulta necesario aplicar diferentes estrategias para reducir el consumo de agua favoreciendo técnicas sostenibles de regadío (García-Tejero y otros [2011a]). El riego deficitario (RD) es una de esas estrategias. Se define como RD la práctica en la que la cantidad de agua aplicada al riego se reduce a una fracción de la evapotranspiración del cultivo (ETc). Según English y Raja [1996], el RD es una estrategia de optimización en la que se permite que los cultivos sufran un cierto grado de déficit hídrico y, consecuentemente, una cierta reducción en la producción. Su aplicación a los cultivos exige un conocimiento exhaustivo de la ETc, de la respuesta del cultivo al déficit hídrico -incluyendo la identificación de los períodos fenológicos críticos-, de los efectos a largo plazo del estrés aplicado durante varias campañas consecutivas y del impacto sobre la producción. Fereres y Soriano [2007] estudiaron el efecto del estrés hídrico (EH) sobre la productividad en los cultivos, comprobando que llega un punto a partir del cual aportes adicionales de riego no se traducen a mejoras en la producción, como se muestra en la figura 4.11. Varios autores han constatado que la respuesta de los cultivos de cítricos al RD presenta una fuerte dependencia con respecto del grado de estrés hídrico soportado por la planta en diferentes estados fenológicos (Goldhammer y Salinas [2000]). Los distintos efectos causados por la estrategia de RD están fuertemente correlacionados con el tiempo de riego (tanto el momento del día como la duración) y el estado fisiológico del cultivo, teniendo en cuenta los períodos de máxima demanda evapotranspirativa, en los que el agua no debe ser escatimada (Chalmers y otros [1986]). Hay diferentes formas de implementar una estrategia de RD, según la forma en la que esa restricción se aplique. El Riego Deficitario Sostenido (RDS) se basa en la aplicación uniforme de restricción hídrica, dependiendo de la demanda 84 Figura 4.11: Relaciones generalizadas entre agua de riego aplicada, EvapoTranspiración y rendimiento del cultivo (adaptada de Fereres y Soriano [2007]). del cultivo. Esta estrategia permite al cultivo desarrollar una adaptación a la situación de estrés. Sin embargo, este enfoque no considera la posible incidencia de periodos críticamente sensibles a la situación de estrés durante el desarrollo fisiológico del cultivo. Fereres y Soriano [2007] establecieron que el estrés hídrico debido a la aplicación de un tratamiento de RDS moderado provoca una reducción en el tamaño de copa y la interceptación de la radiación, lo que afecta negativamente a la fotosíntesis y por lo tanto al crecimiento. El Riego Deficitario Controlado (RDC), por contra, es una estrategia en la que la aplicación de riego depende precisamente de la fisiología del cultivo, y los efectos causados por esta estrategia están fuertemente correlacionados con el tiempo, la duración y el estado fisiológico de la planta, teniendo en cuenta los máximos periodos evapotranspirativos y las fases de crecimiento más críticas en las que el aporte hídrico no debería ser reducido Chalmers y otros [1986]. A pesar de los inconvenientes descritos, el RD ha demostrado en multitud de escenarios y ocasiones ser una política efectiva y sostenible en circunstancias de escasez de agua al minimizar el consumo de agua estabilizando la produc- Modelado, simulación y control del uso del agua en la agricultura 85 ción del cultivo, esto es, maximizando la productividad del agua García-Tejero y otros [2011b]. Concretamente, y dentro de las distintas estrategias de RD, el riego deficitario controlado es la más habitualmente empleada para conseguir ahorro de agua en cultivos leñosos. Aunque los efectos de estas prácticas han sido estudiados ampliamente en distintos cultivos, su aplicación en el sector agrícola todavía se encuentra en sus primeras fases, y quedan por dilucidar ciertas cuestiones: una de ellas consiste en la posibilidad de predecir el efecto que tendrá sobre la producción imponer estrés hídrico en diferentes periodos fenológicos. En Doorenbos y otros [1979] se establece una ecuación para predecir la productividad del cultivo como respuesta al uso de agua, basada a su vez en una ecuación propuesta por Stewart y otros [1977] que definió el factor de proporcionalidad KYentre la pérdida de producción y la reducción de la evapotranspiración asociada. En el ya citado Doorenbos y otros [1979] se proponen varios valores de KY para distintos cultivos para distintas fases de crecimiento del cultivo, así como una KYglobal que abarque todas las fases. De este modo, teniendo en cuenta estas observaciones, se pueden establecer diferentes ecuaciones de modelado para predecir los efectos del déficit hídrico cuando es aplicado en diferentes estados fenológicos. Estos modelos pueden predecir la pérdida de producción (Yi) cuando el suministro de agua no es suficiente para cubrir la demanda total de agua del cultivo. Para ello, parece razonable considerar que la planificación temporal de la aplicación del déficit hídrico durante distintas fases de crecimiento permitiría conseguir una mayor productividad del agua, dado que la respuesta del cultivo depende de la sensibilidad de éste durante ciertos períodos. Jensen [1968] propuso un modelo aditivo basado en transpiración de cultivos para evaluar los efectos del riego en diferentes fases. Este método relacionó la transpiración del cultivo en diferentes etapas con la reducción de la producción; el modelo mostró diferentes relaciones lineales, en las que el valor de la pendiente indicaba el efecto del RD en la producción, permitiendo así comparaciones entre similares niveles de transpiración para diferentes cultivos. En los próximos capítulos se proponen modelos de la reducción de la producción a partir del estrés hídrico, teniendo en cuenta los resultados obtenidos en dos ensayos diferentes: por un lado, un ensayo sobre el efecto a largo plazo (seis años) del RD sobre un cultivo de cítricos; por el otro, un ensayo sobre diferentes tratamientos de RDC a lo largo de tres años consecutivos. De este modo, se presentan funciones de producción, haciendo especial énfasis en las interacciones del estrés hídrico a lo largo del tiempo y en diferentes etapas fenológicas. Asimismo, el efecto sobre cada periodo fenológico (floración, cre- 92 Capítulo 5. Efectos del estrés hídrico a largo plazo Donde τes la constante de tiempo (medida en años), ωla frecuencia angular (rad/año) y los términos bison adimensionales. De este modo, pueden obtenerse diferentes relaciones para cada grado de estrés hídrico. Sin embargo, la expresión general no se puede contemplar sin la variable x. Por lo tanto, se plantea una expresión en variables separables, como se puede ver en 5.4. L(x,t) = X(x)T(t)→L(x0,t) = X(x0)T(t)→L(x,t0) = X(x)T(t0) (5.4) Así, X(x) tiene la expresión polinómica: X(x) = a1x2+a2x+a3(5.5) Que llevado a la ecuación 5.4, queda: L(x,t) = X(t)T(t) = (a1x2+a2x+a3)[b1+b2t+b3e−t/τcos(ωt)] (5.6) Para determinar los valores de los parametros aiybimediante un problema de minimización, es importante fijar ciertas condiciones de contorno. Para empezar, no se produce pérdida de producción cuando el estrés hídrico es nulo en ninguna circunstancia; además, una condición inicial: en cualquier tratamiento, la pérdida en el año 0 (esto es, antes de que empiece la aplicación del estrés), ha de ser nula. Llevadas dichas condiciones a la ecuación 5.6, ésta queda definida tal y como se muestra en la ecuación 5.7: L(x= 0,t),∀t→L(0,t) = a3[b1+b2t+b3e−t/τcos(ωt)],∀t L(x,t = 0),∀x→L(x,0) = (a1x2+a2x+a3)(b1+b3)=0,∀x ⇒       a3= 0 b1=−b3 (5.7) Asimismo, los extremos locales están situados en t=k,k ∈Z⇔ω=π+2nπ rad, n∈N. Por lo tanto, la ecuación 5.6 queda: Modelado, simulación y control del uso del agua en la agricultura 93 L(x,t)=(a1x2+a2x){b2t+b3[e−t/τcos(πt)−1]}(5.8) Dado que esta función fue creada para estimar la pérdida anual de producción, no tiene sentido calcular la función de pérdida L(x,t) en tiempo continuo, por lo que puede ser expresada en forma discreta tal y como muestra la ecuación 5.9. L(x,k) = L(x,t =k)=(a1x2+a2x){b2k+b3[e−k/τ(−1)k−1]},k ∈N(5.9) De este modo, la función de pérdida de producción permite una estimación de la pérdida del año siguiente cuando se ha aplicado el mismo nivel de estrés en los (k-1) años precedentes. Sin embargo, esto excluye la posibilidad de modelar la pérdida de producción en caso de que varíe el nivel de estrés de un año a otro. En un intento de resolver este problema, se propone una transformación en ecuación en diferencias: las ecuaciones se simplificaron bajo la hipótesis de cambios graduales, esto es, x(k)≈X(k−1) ≈X(k−2). Esta aparente contradicción con la idea inicial de variar el nivel de estrés entre varios años se debe a que en la aplicación de los tratamientos teóricos RDS75,RDS65,RDS50 se produjo un error en la aplicación del índice de estrés hídrico, error que fue detectado y medido gracias al uso de los citados pluviómetros analógicos. De este modo, la modificación de la ecuación para utilizar ecuaciones en diferencias responde a la necesidad de adaptarse mejor a los datos del experimento inicial; y dado que estos errores en el IEH son pequeños, la hipótesis de cambios graduales puede aplicarse. De este modo, la ecuación en diferencias queda como se muestra en la ecuación 5.10. L(x,k) = (1−e−1/τ)L(x,k −1)+e−1/τL(x,k −2)+b2(1+e−1/τ)(a1x2+a2x) (5.10) 94 Capítulo 5. Efectos del estrés hídrico a largo plazo Cualquier ecuación en diferencias ha de venir acompañada con ncondiciones de frontera. En este caso, se obtienen dos valores iniciales de la ecuación 5.9. L(x,0) = 0,∀x(5.11) L(x,1) = (a1x2+a2x)[b2−b3(1+ e−1/τ)] (5.12) Para probar nuestro modelo simplificado, ha de compararse con el rendimiento del modelo completo. Esto implica establecer una relación análoga a 5.10 pero considerando la variación en la variable x. Tal y como sugieren J.R. y L.A. [1952], se puede aplicar la transformada Z a la expresión de la ecuación 5.10 con tiempo de muestreo T=1 (donde p=e−1/τ y r.o.c. es la región de convergencia). Z[L(x,t)](z) = L(x,z) = Z[X(x)T(t)](z) = X(x)Z[T(t)](z) L(x,z) = X(x)z−3C1+z−2C2+z−1C1 z−4p2+2z−3p(1 −p)+ z−2(p2+1−4p) +2z−1(p−1)+ 1 r.o.c. :|z|> Max(1,p) (5.13) Donde C1=p2(b2+b3)−pb3;C2=−p2b3+p(2b2+4b3)+b3;C3=−3pb3−b3+b2. Con el uso de la transformada Z, la ecuación 5.13 puede convertirse en una ecuación en diferencias bajo las mismas condiciones iniciales que la ecuación 5.11: L(xk,k) = C1X(k−3)+ C2X(k−2)+ C3X(k−1) +2(1 −p)L(xk−1,k −1) −(p2+1 −4p)L(xk−2,k −2)+ 2p(p−1)L(Xk−3,k −3)−p2L(xk−4,k −4) (5.14) Al llegar a este punto, nos encontramos con dos modelos discretos: el completo 5.14 y el simplificado 5.9. A la hora de decidir cuál de los dos es más Modelado, simulación y control del uso del agua en la agricultura 95 conveniente, en Nash y Sutcliffe [1970] se define el índice de eficiencia de NashSutcliffe para modelos separables. Dicho modelo compara los errores estándar con respecto de la media de dos modelos distintos, ofreciendo información acerca de en qué grado mejora un modelo a otro. La fórmula para calcular ese índice es r2=R2 2−R2 1 1−R2 1= 1 −SSE2 SSE1, donde R2 ies el coeficiente de determinación de mínimos cuadrados R2para los modelos 1 y 2, ySSEies el error estándar con respecto de la media de dichos modelos. 5.2. Resultados y discusión 5.2.1. Relaciones hídricas y respuesta fisiológica del cultivo La tabla 5.1 muestra las condiciones meteorológicas y las dotaciones hídricas recibidas por cada uno de los tratamientos estudiados durante seis campañas de riego (2005-2010), con valores diarios de la evapotranspiración del cultivo de 5,34, 4,54, 3,20, 3,78, 3,17 y 3,26. Las precipitaciones se distribuyeron fundamentalmente entre el final del otoño y el principio de la primavera, mientras que en el periodo de máxima demanda evapotranspirativa fueron prácticamente despreciables. Las cantidades de riego en los tratamientos estresados se aproximaron a funciones de la evapotranspiración del cultivo, tal y como estaba previsto en su diseño, y el ahorro medio en agua de riego con respecto al tratamiento control para los tratamientos RDS 50, 65 y 75 fue del 46%, 34% y 23% respectivamente. Se pudo observar que el aporte de agua (incluyendo tanto riego como precipitaciones) en el tratamiento control C100 fue superior a la evapotranspiración del cultivo, aunque este hecho fue el resultado de la aparición de tormentas tanto al principio como al final de las campañas de riego, coincidiendo con el principio de la primavera y el final del otoño. 96 Capítulo 5. Efectos del estrés hídrico a largo plazo Tabla 5.1: Agua aplicada, Balance hídrico y tratamientos de riego. El periodo de riego se mide en días del año (DDA). Año 2005 2006 2007 2008 2009 2010 Periodo de riego 150-280 110-275 110-309 116-284 110-329 115-277 ETc(mm) 700 754 641 639 698,7 547,8 Precip. (mm) 28 119 115 117 67,2 89,4 RDS50 Riego (mm) 372 348 356 345 356,3 257,5 EH 0,53 0,45 0,56 0,54 0,51 0,47 RDS65 Riego (mm) 465 435 445 409 444,4 345 EH 0,66 0,56 0,69 0,64 0,64 0,63 RDS75 Riego (mm) 533 501 512 472,9 517 399,9 EH 0,76 0,67 0,80 0,74 0,74 0,73 C100 Riego (mm) 696 653 668 632.6 680.1 531.4 EH 0,99 0,88 1,04 0,99 0,97 0,97 EH: Estrés hídrico; RDS: Riego deficitario sostenido. 5.2.2. Respuesta productiva del cultivo a riego deficitario sostenido y funciones de productividad del agua anuales. Los efectos del riego deficitario sobre la producción frutícola a lo largo de seis años pueden observarse en la tabla 5.2 Tabla 5.2: Producción media de frutos (kg/árbol) a lo largo del experimento. Diferentes letras muestran diferencias significativas entre tratamientos (p < 0,05). Tratamiento Año 1 Año 2 Año 3 Año 4 Año 5 Año 6 RDS50 88,3a119,6a76,0a120,5a47,1a75,7a RDS65 100,7a126,7a81,7a121,0a48,0a78,0ab RDS75 99,1a113,9a82,0a121,0a60,6b88,9ab C100 111,3a124,0a95,7a132,5b68,8b99,8b RDS: Riego deficitario sostenido; C: Tratamiento control Durante los primeros años del experimento, entre 2005 y 2008, no se encontraron diferencias significativas entre el tratamiento control y los tratamientos deficitarios. Es en el periodo siguiente (2008-2010) cuando empezaron a obser- Modelado, simulación y control del uso del agua en la agricultura 97 varse diferencias significativas en la producción, de forma especialmente apreciable en los tratamientos de mayor estrés, RDS50 y RDS65. Finalmente, entre 2009 y 2010, no se apreciaron diferencias significativas entre el tratamiento control y el RDS75. La pérdida de la producción de un árbol, en referencia a la producción del tratamiento control cuando ese árbol está sometido a los distintos tratamientos de déficit hídrico varió entre el 16% para los árboles sometidos al tratamiento RDS50 y el 9% para el tratamiento RDS75. La figura 5.3 muestra la relación entre pérdida de producción (1−Y ieldi Y ieldC100 ) y el estrés hídrico aplicado (1−ETci/ETcC100) año a año. Según esos resultados, en el año 2005 no se produjo un efecto significativo sobre la producción, aunque existió una relación directa entre ambas magnitudes. Figura 5.3: Coeficientes Kyde funciones de productividad del agua (CWPF) lineales para cada año relacionando la pérdida de producción con la tasa de evapotranspiración aplicada (ETc). 98 Capítulo 5. Efectos del estrés hídrico a largo plazo Esta relación puede verse en la tabla 5.3 donde se muestran las funciones de productividad del agua CWPF para ese año con Ky= 0,39 y un coeficiente de determinación R2= 0,90. Tabla 5.3: Valores estimados de Kyy coeficiente de determinación R2para funciones CWPF en función de ETcyψstem. CWPF Año 1 Año 2 Año 3 Año 4 Año 5 Año 6 1−Yi YC100 =Ky ETci ETcC100 Ky0,39 0,01 0,49 0,24 0,68 0,48 R20,90∗∗ ns 0,87∗∗ 0,84∗∗ 0,89∗∗ 0,94∗∗ **: significativo con p < 0,01; *: significativo con p < 0,05 ns: no significativo En el año siguiente, 2006, la relación fue despreciable, con una función CWPF prácticamente plana, pendiente cercana a cero, sin diferencias importantes entre los tratamientos. Al llegar a 2007 las relaciones volvieron a ser significativas, con Ky= 0,39 y R2= 0,87; sin embargo, el efecto sobre la producción volvió a ser despreciable. Sólo a partir del siguiente año, 2008, las diferencias entre los distintos tratamientos RDS empiezan a ser significativas, y esa significación sigue creciendo en los tres años siguientes hasta el fin del estudio en 2010. Los resultados del experimento mostraron una clara tendencia en la respuesta del cultivo al estrés hídrico a lo largo de varios años consecutivos, sugiriendo que los efectos no son similares para diferentes estaciones y que la Kyno muestra una tendencia al crecimiento a lo largo del tiempo. Es más, existe una oscilación en los valores de Ky, con años en los que los valores de esta magnitud son elevados (2005, 2007 y 2009, los años impares) y años con Kymás bajos (los años pares correspondientes). Es por eso que, a la hora de modelar los efectos del riego deficitario sobre la producción en cítricos, se ha de tener en cuenta dos parámetros: el índice de estrés hídrico y la duración de dicho estrés. Varios autores que han estudiado funciones CWPF en cultivos anuales han mostrado relaciones lineales entre las pérdidas de producción y el agua aplicada (Stewart y Hagan [1973]; Hanks [1974]; Hanks [1983]; Hexem y Heady [1978]). Estas relaciones pueden resultar útiles para el diseño de estrategias de riego, eligiendo los tiempos de riego y la cantidad de agua aplicada, permitiendo además comparar eficiencias relativas del uso del agua (Al-Jamal y otros [2000]). Sin embargo, resulta complicado establecer una única función CWPF para un cultivo leñoso perenne, ya que en ese caso entran en juego diversos factores que Modelado, simulación y control del uso del agua en la agricultura 99 provocan variaciones en los valores de Ky. Entre esos factores destacan, como se ha visto, la variabilidad temporal, pero también la espacial (depende de la parcela) y los estados fenológicos del cultivo en los que se aplica el déficit hídrico. En este sentido, Clumpner y K. [1987] y Sammis [1981] obtuvieron resultados similares a los de este estudio. De acuerdo con esta argumentación, para cultivos anuales, Kyvariará según el lugar pero los valores deberían ser similares año a año, viéndose afectadas tan sólo por las condiciones meteorológicas (Al-Jamal y otros [2000]). De este modo, si se acepta la hipótesis de que las condiciones climáticas fueron las mismas durante todo el periodo del ensayo, la variación interanual de Kypodría estar relacionada con el efecto acumulativo del estrés hídricos, y la capacidad de respuesta de un árbol sometido a dicho estrés acumulado a lo largo de varios años. 5.2.3. Funciones de pérdida de producción De acuerdo con los resultados anteriores, puede definirse una función de pérdida de producción en frutos L(x,t), donde xes el nivel de estrés hídrico y t es el tiempo. Esta función permitiría una estimación de la pérdida de producción en un año kcuando se aplique un cierto nivel de estrés xen ese año y en los k−1 años anteriores. Aplicando el cálculo de la minimización de la función definida en la ecuación 5.2, se obtiene: L(x,t) = (0,4127x−0,1553x2){0,1871t−0,4731[e−t/12,11cos(πt)−1]}(5.15) La figura 5.4 muestra la curva tridimensional que sigue la pérdida de producción como función de los años de observación y el estrés hídricos. Para la evaluación de esta función de pérdida de producción, se agruparon los datos según su nivel de etrés, analizando la bondad del ajuste en tres escenarios diferentes de riego deficitario sostenido. L(x=xi,t) = k(xi){b2t+b3e−t/τ[cos(πt)−1]}(5.16) 100 Capítulo 5. Efectos del estrés hídrico a largo plazo Donde xison los tres escenarios diferentes estudiados. Figura 5.4: Curva de superficie para la función de pérdida de producción y=L(x,t); comparación del modelo con los resultados experimentales. RMSE: Raíz del error cuadrático; SSE: Suma de los residuos al cuadrado; r2: Coeficiente de determinación. La tabla 5.4 muestra las relaciones de significación calculadas para los tres niveles de estrés hídrico considerados para el ajuste en tiempo continuo, considerando esas relaciones definidas por el estrés hídrico durante seis años consecutivos. Se puede observar que dichas relaciones vendrán definidas por tres funciones oscilantes en las que los efectos del estrés no son constantes a lo largo del tiempo del ensayo. Sin embargo, se pudo observar cambios en dichos efectos en forma de recuperación parcial seguidos de años con pérdidas significativamente mayores en la producción (figura 5.4). También puede observarse que la amplitud de esas oscilaciones va decreciendo progresivamente, aunque ese efecto de atenuación es mayor para niveles de estrés menores. El modelo propuesto define una serie de funciones para calcular la función de pérdida de producción en función del tiempo, asumiendo que éste es una variable continua. Sin embargo, los efectos del estrés hídrico sobre la producción en un cultivo frutal como el de los cítricos sólo pueden ser evaluados al Modelado, simulación y control del uso del agua en la agricultura 101 Tabla 5.4: Cálculo de las estadísticas de la bondad del ajuste en los escenarios de estrés hídrico analizados. Índice de estrés hídrico RMSE SSE r2 (kg/kg) (10−3[Kg/Kg]2) Tiempo continuo RDS50 0,039 1,509 0,88∗∗ RDS65 0,065 4,280 0,66∗∗ RDS75 0,018 0,323 0,85∗∗ Tiempo discreto simplificado RDS50 0,037 1,40 0,89∗∗ RDS65 0,064 4,1 0,68∗∗ RDS75 0,018 0,315 0,85∗∗ Tiempo discreto completo RDS50 0,023 0,513 0,96∗∗ RDS65 0,046 2,144 0,83∗∗ RDS75 0,031 0,978 0,54∗∗ **: significativo con p < 0,01; RMSE: Raíz del error cuadrático. SSE: Suma de los resíduos al cuadrado. r2: Coeficiente de determinación. final de la temporada, por lo que la variable tiempo debería ser considerada como discreta; sólo año a año se podrá decidir qué estrategia de gestión del estrés hídrico se ha de seguir y observar sus efectos sobre la producción. A la vista de este hecho, no tendría sentido seguir calculando la función L(x,t) en tiempo continuo. Se definió un modelo en tiempo discreto para estimar la pérdida de producción anual cuando se aplicaba el mismo nivel de estrés (o un nivel similar) en los años anteriores. Tras llevar a cabo el proceso de minimización del error, se obtuvieron los siguientes resultados: a1=−0,1567 a2= 0,4342 b2= 0,1669 b3=−0,4935 τ= 12,1214 108 Capítulo 6. Efectos del estrés hídrico a medio plazo sudeste. El área sombreada del suelo y la zona de riego forman el 30% y el 17% de la superficie total. El suelo de la zona del ensayo es Fluvisol típico FAO98 [1998], con textura franco-arcillo-arenosa (350 g/kg de arena; 400 g/kg de limo; 250 g/kg de arcilla) y contenido de materia orgánica por debajo de los 15 g/kg. Tiene una profundidad de 1,5m, con capacidad de campo y punto de marchitez permanente de 235 mm/m y 100 mm/m respectivamente. La climatología local es mediterránea, con una ETomedia de 1500 mm/año y precipitaciones anuales de 475 mm, distribuidas fundamentalmente entre finales del otoño y principios de la primavera, con una gran variabilidad interanual. Las temperaturas se ven sometidas a grandes variaciones a lo largo del año, con temperaturas medias que van desde los 10◦Cdel invierno a los 35◦Cdel verano, superando a menudo los 40◦Cdurante las horas de mayor efecto de la radiación solar. 6.1.2. Diseño experimental. Tratamientos RDC El diseño experimental fue de bloques al azar, con cinco repeticiones por tratamiento. La unidad experimental (360 m2) consistía tres filas con cinco árboles por fila, y los tres árboles centrales (en lo sucesivo, fila central) fueron los escogidos para medir la productividad del cultivo frutal. La distribución puede verse en la figura 6.1. Tabla 6.1: Definición de los tratamientos de riego. Tratamiento Floración Crecimiento Maduración RDI10,55 0,70 0,55 RDI20,55 0,70 0,70 RDI30,55 0,55 0,70 RDI40,70 0,70 0,55 C100 Regado permanentemente al 100% de la ETc Los tratamientos de riego se establecieron en función de diferentes índices de estrés, definidos a través de la relación entre el volumen de agua aplicado en cada tratamiento y el valor de la evapotranspiración real del cultivo (ETc), teniendo en cuenta la etapa fenológica del cultivo (floración, crecimiento o maduración), según se relacionan en la tabla 6.1. Se instalaron contadores volumétricos para registrar el volumen de agua aplicada a cada tratamiento. Modelado, simulación y control del uso del agua en la agricultura 109 Figura 6.1: Distribución de los tratamientos de riego en el ensayo de RDC. Las distintas fase consideradas se corresponden con los principales periodos fenológicos en el cultivo del cítrico. La primera fase abarca la floración y el establecimiento del fruto; la segunda fase, los meses de julio y agosto en los que la fruta alcanza el tamaño máximo; finalmente, durante la fase de maduración, se producen diversos cambios químicos en el fruto, que provocan un cambio de color y un incremento en el total de sólidos solubles. Los valores estacionales de la ETcfueron determinados empleando la ecuación de Doorenbos y Pruitt (Doorenbos y Pruitt [1977], ecuación 5.1), recogiendo datos meteorológicos en una estación automatizada situada en las proximidades del cultivo. Los coeficientes del cultivo variaron entre 0,50 entre marzo y mayo, 0,55 entre junio y octubre y de nuevo 0,50 en noviembre y diciembre. Al final de cada estación se analizaron para su cálculo los valores de la producción de fruto, la productividad del agua y la productividad del agua de riego. Se llevó a cabo un análisis exploratorio y descriptivo de la producción seguido de un análisis con un análisis de separación de media usando el test de significación estadística honesta (HSD) de Tukey Linton y Harder [2007]. 110 Capítulo 6. Efectos del estrés hídrico a medio plazo El objetivo principal fue determinar una expresión para determinar la pérdida de producción en función del nivel de estrés hídrico aplicado. La función óptima se determinó mediante la minimización del error cuadrático medio (RMSE), tal y como se mostró en la ecuación 5.2 en el capítulo 5. Donde ˆ f(x) es cada una de las posibles funciones de producción, x=x(t) es el nivel de estrés hídrico aplicado (que puede ser dependiente del tiempo) e yi es la pérdida de producción medida tras la aplicación del nivel de estrés xi. 6.1.3. Modelando la función de producción del cultivo Se determinaron diferentes relaciones entre el estrés hídrico aplicado y la pérdida de producción, teniendo en cuenta que los diferentes niveles de estrés hídrico fueron definidos en diferentes fases: floración, crecimiento de fruto y maduración. Si se parte de la hipótesis de que el estrés hídrico en una fase no tiene efectos sobre las siguientes: f(x) = f(xf,xg,xr) (6.1) Donde f(x) es la función de pérdida de producción y los índices de estrés en las fases de floración, crecimiento de fruto (growing) y maduración (ripening) vienen dados por las variables xf,xg,xrrespectivamente. Se propusieron funciones cuasi-lineales, cuadráticas y cúbicas, obteniendo la máxima precisión con las últimas. Para resolver el problema de minimización se emplearon las funciones de minimización del software matemático MATLAB®(The Mathworks Inc., Natick, Ma., EE.UU.) y el algoritmo de programación secuencial cuadrático Fletcher y Powell [1963] Gill y otros [1981], que nos proporcionaron las funciones más adecuadas. La ecuación general para cada fase se definió como se muestra en la ecuación 6.2. fk(xk) = akx3 k+bkx2 k+ckxk(6.2) Modelado, simulación y control del uso del agua en la agricultura 111 Considerando que cuando no se produce estrés hídrico, no hay pérdida de la producción, es razonable asumir que la función fk(xk) no tenga ordenada en el origen: fk(0) = 0. 6.1.4. Modelo aditivo La ecuación general que incluye las tres fases es la suma de las tres ecuaciones previas definidas en la ecuación 6.2. Llevadas a la ecuación 6.1, se obtiene: f(x) = f(xf,xg,xr) = ff(xf)+fg(xg) +fm(xm) (6.3) Se tuvieron en cuenta cuatro restricciones para cada función fk(xk): 1. En el rango de interés, xk∈[0,1], la función debe ser monótona creciente: f′ k(xk)≥0,∀xk∈[0,1] →3akx2 k+2bkxk+ck≥0 (6.4) 2. La productividad no puede reducirse en un factor mayor que 1: fk(xk)≤1,∀xk∈[0,1] (6.5) 3. Por definición, la función vale 0 en el origen fk(0) = 0 (6.6) 4. Para que el modelo aditivo sea coherente, la suma de las tres componentes no puede ser mayor que uno en ningún momento (esta condición hace la condición 6.5 redundante). Como, por la condición 6.4 las funciones son monótonas crecientes: ff(1)+ fg(1)+ fm(1) ≤1→f(1,1,1) = ∑ k=f ,g,r ak+bk+ck≤1 (6.7) 112 Capítulo 6. Efectos del estrés hídrico a medio plazo Dadas esas cuatro restricciones, la función a minimizar es: Min :RMSE =Er·ET r Er= (AX +−B)⊘B A= [kn ij],B = [yi]T;i= 1..,4;j=f ,g,r;n= 1..,3 X=[afagarbfbgbrcfcgcr] (6.8) Lo que forma un problema de minimización en forma matricial donde la incógnita es el vector X. kij son los niveles de estrés aplicados en la fase j del experimento i; la expresión M⊘Nimplica el cociente de Hadamard (elemento a elemento) de dos vectores M y N Byung-Gook y Yunbeom [1998]. Finalmente, y siguiendo a Moriasi y otros [2007], se empleó el índice de eficiencia de NashSutcliffeR2para cuantificar la eficiencia del modelo resultante Nash y Sutcliffe [1970]. 6.1.5. Modelo por ajuste con términos cruzados La principal hipótesis de partida de este modelo es que el estrés inducido en dos periodos fenológicos diferentes debería tener un efecto mayor que su aplicación en tan sólo una fase, pero podría tener un efecto más débil que la simple suma de dos componentes separados. Dado que fes una función de pérdida de producción, si el estrés se induce en la primera fase, la producción de fruto se reducirá; consecuentemente, el estrés hídrico aplicado en periodos posteriores no podrá tener efecto sobre la parte de la producción que ya se ha perdido. f(xf,0,0) = ff(xf)→P(t1) = 1−f(xf,0,0) = 1 −f(xf) f(xf,xg,0) = ff(xf)+fg(xg)P(t1)→P(t2)=1−f(xf,xg,0) f(xf,xg,xr) = ff(xf) +fg(xg)P(t1)+ fr(xr)P(t2) (6.9) Modelado, simulación y control del uso del agua en la agricultura 113 Donde P(t) es una función que describe la producción resultante de aplicar estrés hídrico (es decir, el complementario de las pérdidas); t1,t2son los instantes de tiempo al final de las fases de floración y crecimiento respectivamente. Por lo tanto, P(t1),P (t2) son las producciones resultantes de aplicar estrés hídrico en las dos primeras fases. Obviamente, si tomamos t0como el instante de tiempo anterior a la fase de maduración, P(t0) = 1. La ecuación 6.9 se puede desarrollar: f(xf,xg,0) = ff(xf)+ fg(xg)P(t1) = ff(xf)+ fg(xg)−ff(xf)fg(xg) f(x) = f(xf,xg,xr) = ff(xf)+fg(xg)P(t1) +fr(xr)P(t2)→ f(x) = ff(xf) +fg(xg)+ fr(xr)−ff(xf)fg(xg)−ff(xf)fr(xr) −fg(xg)fr(xr)+ ff(xf)fg(xg)fr(xr) (6.10) Esta ecuación puede escribirse en forma resumida tal y como muestra 6.11. f(xf,xg,xr) = ∑ i=f ,g,r fi(xi)+ ∑ i=f ,g,r;i,j fij(xi,xj) +ff grx(xf,xg,xr) (6.11) Donde las fk(xk) pueden calcularse con las restricciones previamente propuestas, y las fij(xi,xj) y ff gr(x) se pueden calcular según la ecuación 6.12. fij(xi,xj) = −fi(xi)fj(xj) ff gr(x) = ff(xf)fg(xg)fr(xr)(6.12) Si se garantiza el cumplimiento de esas restricciones, el efecto del estrés hídrico en dos fases ha de ser mayor (o, en el caso degenerado, igual) que el efecto del estrés en una cualquiera de ellas, lo que viene garantizado por la restricción 6.5: 0 ≤fi(xi)≤1,∀xi∈[0,1]. 114 Capítulo 6. Efectos del estrés hídrico a medio plazo Además, la superposición de dos de estos efectos ha de dar como resultado una función monótona creciente, lo que a su vez viene garantizado por las condiciones 6.4 y 6.5. fij +fj+fi=fi+fj−fifj≤1,∀(xi,xj)∈[0,1]2→fj(1−fi)≤1−fi 1−fi≥0,∀xi∈[0,1] →fj(xj)≤1,∀xj∈[0,1] (6.13) Sumando los términos cruzados de las tres fases (floración, crecimiento de fruto y maduración), el efecto del estrés hídrico ha de ser mayor que el efecto del estrés aplicado tan sólo sobre dos de dichas fases cualesquiera. Esto nos lleva a la ecuación 6.14, garantizada por la restricción 6.5. fi+fj+fk−fifj−fifk−fjfk+fifjfk≥f i +fj−fifj→fk−fifk−fjfk+fifjfk≥0 fk(xk)≥0,∀xk∈[0,1] →1−fi−fj+fifj≥0→fi+fj−fifj≤1,∀(xi,xj)∈[0,1]2 (6.14) Como indicaba la condición 6.5, la pérdida de producción tiene una cota superior de 1, como muestra la ecuación 6.15. f(x) = ff+fg+fr−fffg−fffr−fgfr+fffgfr≤1,∀(xf,xg,xr)∈[0,1]3 1−fg(xg)≥0,∀xg∈[0,1] →ff(1−fg)+fr(1 −fg)+fffr(1 −fg)≤1−fg→ ff+fr+fffr≤1,∀(ff,fr)∈[0,1]2 (6.15) Desgraciadamente, no es posible transformar todas las restricciones del problema en condiciones lineales. La condición 6.4 (dfk(xk) dxk≥0,∀xk∈[0,1]) ha de ser satisfecha en todos los puntos del intervalo [0,1], dando lugar a infinitas restricciones lineales. Para evitar este problema, se escogieron una serie de puntos representativos (por ejemplo, 0,0,5,1) lo bastante numerosos para evaluar la restricción en dichos puntos. Las restricciones sobre esos puntos se almacenaron Modelado, simulación y control del uso del agua en la agricultura 115 en un archivo de script Matlab®(.m), dando lugar a un total de 20 restricciones para las nueve variables del problema. El problema de minimización propuesto para resolver este modelo sigue estando definido por la ecuación 6.8, pero las matrices han cambiado tanto que el problema debe ser redefinido: Min :RMSE =Er·ET r Er=ˆ F⊘B B=[y1y2y3y4]T Fi=[aik3 ij bik2 ij cikij];j= 1..,4;i=f ,g,r P1=[1 1 1 1]T−Ff P2=[1 1 1 1]T−(Ff+Fg⊙P1) ˆ F=Ff+Fg⊙P1+Fr⊙P2 (6.16) Donde M⊙Nimplica el producto de Hadamard (término a término) de dos vectores Byung-Gook y Yunbeom [1998]. Finalmente, se evaluó el coeficiente de determinación R2. Además, para comparar los dos modelos (el aditivo y el de términos cruzados) se emplea el índice de mejora r2que permite evaluar el grado de mejora que aporta un modelo sobre otro (Nash y Sutcliffe [1970]), actuando como descripción cuantitativa de hasta qué punto mejora la fiabilidad del modelo añadir los términos cruzados. 116 Capítulo 6. Efectos del estrés hídrico a medio plazo 6.2. Resultados y discusión 6.2.1. Producción y tratamientos de RDC El efecto del déficit hídrico sobre la producción de fruto y la productividad del agua fueron evidentes, con ciertas diferencias significativas entre tratamientos que se pueden observar en la tabla 6.2. Tabla 6.2: Efectos del RDC sobre la producción y la productividad del agua en los tres periodos de monitorización de un año. Trats. Riego Riego+Lluvia Producción Reducción PA PAR (mm) (mm) kg/árbol (%) (g/L) (g/L) Año 1 RDC1180 513 125,2bc 18,7 8,1c23,2b RDC2209 542 130,9b14,9 9,7a25,1a RDC3174 507 121,9c20,8 8,0c23,3b RDC4215 548 142,2a7,6 8,6b22,0b C100 311 644 153,9a8,0c16,5c Año 2 RDC1207 470 154,5b15,3 10,9b24,9b RDC2212 475 153,2b16,0 10,8b24,1b RDC3183 446 147,7c19,0 11,0b26,9a RDC4215 478 169,6ab 7,0 11,8a26,3a C100 318 581 182,4a10,4c19,1c Año 3 RDC1214 539 129,8b17,8 8,0b20,2b RDC2219 544 121,4b23,1 7,2c18,4c RDC3185 510 117,2c25,8 7,7c21,1a RDC4219 544 141,3a10,5 8,6a21,5a C100 322 647 157,9a8,1b16,3d Las diferentes letras (a,b,c) indican diferencias significativas con p < 0,05 según el test HSD de Tukey PA: Productividad del Agua; PAR: Productividad del Agua de Riego Ciertos tratamientos en los que se aplicaron cantidades similares de agua de riego proporcionaron distintas respuestas productivas, lo que puede correlacionarse con la estrategia de riego seguida para cada tratamiento; tal y como se Modelado, simulación y control del uso del agua en la agricultura 117 evidencia en García-Tejero y otros [2008], mostrando la importancia del periodo fenológico en que se aplica el estrés hídrico. En términos productivos, los mejores resultados se consiguieron con el tratamiento RDC4, en el que se aplicó un severo estrés hídrico durante la fase de maduración: el cultivo recibió durante ese periodo un 55% de la ETc; en cambio, durante los otros dos periodos las restricciones se limitaron a proporcionar al cultivo un 70% de la ETc. Otros tratamientos de RDC con similares cantidades de agua aplicada mostraron respuestas productivas positivas, pero en ningún caso tan positivas como en dicho tratamiento. En este sentido, la cantidad de agua aportada al cultivo en el segundo tratamiento RDC2, en el que se regó hasta lograr el 55% de la ETcdurante la fase de floración y el 70% en las fases de crecimiento de fruto y su maduración, fue similar al del cuarto tratamiento, pero la producción se vio considerablemente más restringida (llegando a perder un 18% de media sobre la producción del tratamiento control). Estos resultados resaltan la importancia del periodo fenológico en el que se aplica el estrés hídrico. El tercer tratamiento RDC3proporcionó un ahorro de agua de aproximadamente 130 mm y una reducción de la producción del 22% con respecto al control, mientras que el RDC1, con una media de 117 mm de ahorro de agua, alcanzó valores relativos de producción similares a los obtenidos con el RDC2. Estos tratamientos se vieron sometidos a las mismas restricciones de riego durante los periodos fenológicos de floración y crecimiento respectivamente, pero el tratamiento RDC2tuvo una recuperación parcial durante la fase de maduración; esta recuperación parcial, contrariamente a lo que cabría esperar, no supuso un incremento significativo en la producción con respecto al RDC1. Por lo tanto, el esfuerzo de riego durante la fase de maduración no puede ser considerado como especialmente significativo en comparación con las otras dos fases. En término medio el tratamiento RDC4mostró una reducción del 8%, con ahorros de agua próximos a 101 mm. Es digno de reseñarse que las restricciones más severas de este tratamiento se llevaron a cabo durante la fase de maduración, lo que resulta coherente con la hipótesis anterior y resalta a su vez la importancia de las fases de floración y crecimiento en cuanto a la producción final. El efecto global de los efectos del estrés hídrico sobre la producción están fuertemente correlacionados con los periodos fenológicos en los que se aplicó el estrés, según los resultados de González-Altozano y Castel [2000]. Por otro lado, los mismos autores en González-Altozano y Castel [2003], estudiando cí- 124 Capítulo 6. Efectos del estrés hídrico a medio plazo La tabla 6.3 muestra los resultados de repetir estos cálculos con funciones polinómicas de grado 2 y grado 4. Tabla 6.3: Errores e índices de eficiencia para modelos con funciones polinómicas de distintos grados. Grado de las funciones 2 3 4 Sin términos cruzados ERMS 0,0597 0,0304 0,0171 R20,7922 0,9959 0,9988 Con términos cruzados ERMS 0,0105 0,0105 0,0029 R20,8098 0,9995 1 r20,0847 0,8700 1 ERMS: raíz del error cuadrático medio R2: coeficiente de eficiencia de Nash-Sutcliffe r2: índice de mejora del modelo con términos cruzados A la vista de los resultados de la tabla: 1. Las funciones de grado 2 proporcionaron unos resultados con bajo índice de eficiencia (R2= 0,80). Los términos cruzados sólo tuvieron un efecto marginal a la hora de mejorar esa eficiencia (r2= 0,0847). 2. Las funciones de grado 3 redujeron a la mitad el error RMS y proporcionaron un índice de eficiencia próximo a la unidad. La adición de los términos cruzados mejoró aún más esa eficiencia (r2= 0,87). 3. Las funciones de grado 4 redujeron el error RMS a la mitad una vez más llegando a una eficiencia igual a la unidad al producirse el añadido de los términos cruzados. El grado 4 sería óptimo pero el grado 3 es, incluso desde el punto de vista más exigente, suficientemente bueno y reduce las posibilidades de un exceso de condicionamiento del problema. Por eso se escoge el grado 3. La tabla 6.4 muestra los resultados esperados según el modelo con términos cruzados con funciones de grado 3 y los resultados obtenidos en el ensayo. Se puede observar que las diferencias entre ambos son pequeñas, entre −0,05 y 0,04 de reducción de producción. Estos resultados, que pueden observarse gráficamente en la figura 6.5, sugieren que el modelo con funciones de grado 3 tiene suficiente precisión como Modelado, simulación y control del uso del agua en la agricultura 125 Tabla 6.4: Tasas de pérdida de producción (kg/kg) medidas (med.) y previstas (prev.) con respecto a la producción del tratamiento control según tratamiento y año. Tratam. Campaña 1 Campaña 2 Campaña 3 Med. Prev. Error Med. Prev. Error Med. Prev. Error RDC10,19 0,21 0,02 0,15 0,19 0,04 0,18 0,20 0,02 RDC20,15 0,19 0,04 0,16 0,17 0,01 0,23 0,18 −0,05 RDC30,21 0,23 0,02 0,19 0,21 0,02 0,26 0,22 −0,04 RDC40,08 0,11 0,03 0,07 0,11 0,04 0,11 0,12 0,01 RDC40,21 0,23 0,02 0,19 0,21 0,02 0,26 0,22 −0,04 Figura 6.5: Comparación gráfica entre los resultados del modelo (mod) y los experimentales (exp) para los cuatro tratamientos a lo largo de tres campañas. para estimar una pérdida de producción cuando el estrés hídrico es aplicado durante dos o más periodos fenológicos. 126 Capítulo 6. Efectos del estrés hídrico a medio plazo Se han definido diferentes factores de respuesta del cultivo Kysegún la metodología propuesta por varios autores: Stewart y Hagan [1973], Azaiez [2008], Houndari y Azaiez [2001]. Sin embargo, Stewart y Hagan [1973] mostraron diferentes relaciones entre producción y evapotranspiración del cultivo, desarrollando una ecuación multiplicativa para la producción del cultivo en función del agua de riego aplicada. Por otro lado, Hexem y Heady [1978] definieron una ecuación general para definir las funciones de producción del agua como polinomios de tercer orden. Kipkorir y otros [2002] encontraron un ajuste significativo con polinomios de segundo y tercer orden para cultivos de cebollas y maíz respectivamente. En este sentido, en Severini y Cortignani [2008] se propone un modelo de programación matemática positiva para estudiar la viabilidad de las estrategias de riego deficitario, exponiendo que esas técnicas describen de forma eficaz el impacto negativo de la limitación de la disponibilidad de agua en los resultados económicos de explotaciones agrícolas. Además, en Azaiez [2008] se desarrolla un modelo integral de programación dinámica para analizar los efectos del riego deficitario para incrementar el área regada a costa de reducir la producción por unidad de superficie, considerando que los efectos del riego deficitario están fuertemente correlacionados con la cantidad de agua aplicada para cada periodo fenológico. 6.3. Conclusiones Cuando se aplica déficit hídrico en una fase determinada, la respuesta del cultivo depende tanto de la magnitud de ese estrés como de su duración. Estos efectos son fuertemente dependientes del cultivo, tanto en especie como en variedad: la sensibilidad del cultivo al estrés hídrico varía, incluso la importancia de cada periodo fenológico -que es el principal objeto de este ensayo-. Generalmente, los efectos de la escasez de agua en diferentes periodos interaccionan entre sí de forma compleja, aunque en la mayoría de los casos se pueden introducir simplificaciones mediante la hipótesis de que los efectos del estrés son independientes para cada periodo. En consecuencia, los efectos combinados del estrés hídrico a lo largo de varios periodos se tienden a evaluar considerando esos efectos como aditivos o multiplicativos; dichos modelos re- Modelado, simulación y control del uso del agua en la agricultura 127 quieren calibración local, como se indica en Rao y otros [1988], dado que esas funciones de producción son específicas del ensayo. En un contexto de escasez de agua, estrategias de riego deficitario pueden resultar clave en posibles futuros escenarios de escasez de agua, dado que podrían permitir la supervivencia de agro-ecosistemas en condiciones limitadas. Funciones como las obtenidas en el presente capítulo mejoran nuestro conocimiento respecto a la relación entre el estrés hídrico y la producción de cítricos, y por lo tanto sobre sus implicaciones a la hora de buscar una gestión sostenible de los recursos hídricos disponibles. Según los resultados del estudio, el nuevo modelo propuesto para explicar el efecto combinado del riego deficitario en varios periodos sería un modelo híbrido aditivo-multiplicativo, considerando que el efecto del estrés hídrico en periodos anteriores determina el efecto que tendrá el mismo estrés sobre periodos posteriores. Nuestros modelos indican que la función de productividad del agua para el cultivo sometido a estrategias de riego deficitario tendría una relación cuasi-lineal para déficits hídricos por debajo del 40% de la ETc. Concluimos que los modelos multiplicativos considerando la adición de términos cruzados son aplicables en un rango amplio de condiciones de estrés, especialmente cuando se aplica estrés hídrico en dos o más periodos fenológicos. Simultáneamente, los modelos simples sin considerar esos términos cruzados pueden ser válidos cuando sólo se aplica el estrés en uno de los periodos. Asimismo, analizando la relación entre el estrés hídrico y la pérdida de producción, deducimos que los periodos de floración y crecimiento de fruto (especialmente el primero de ellos) son más sensibles al riego deficitario que el periodo de maduración. De este modo, los modelos desarrollados pueden ser considerados como un paso previo para conseguir una gestión más sostenible y precisa de los recursos hídricos bajo estrategias de riego deficitario en cítricos, especialmente en situaciones de escasez de agua, que son cada vez más frecuentes en zonas áridas y semi-áridas como el sur de España. Capítulo 7 Planificación de una estrategia de riego deficitario en cítricos a medio y largo plazo 7.1. Introducción Como se ha visto, es posible crear una serie de modelos matemáticos para prever la respuesta productiva del cultivo a diferentes grados de estrés hídrico gracias a las funciones de productividad del agua en el cultivo o CWPF. El mero hecho de mejorar nuestro conocimiento acerca de la fisiología del cultivo ya es más que suficiente justificación para emprender la tarea de elaborar esos modelos, pero eso no significa que debamos detenernos ahí. El empleo de técnicas de control y de optimización nos permitiría optimizar los recursos hídricos disponibles, maximizando así la producción y los beneficios. En primer lugar, es preciso afrontar el primer y principal problema de un agricultor que desarrolle su actividad en una zona de clima árido o semi-árido en lo referente al riego: las restricciones. De acuerdo con las observaciones meteorológicas disponibles, se puede esperar significativos incrementos en las 129 130 Capítulo 7. Planificación de estrategias de riego temperaturas globales y descensos en las precipitaciones, lo que provocará de forma inmediata una reducción de la disponibilidad de recursos hídricos (Turral y otros [2011]). Dado que el 70% del agua dulce disponible del mundo se consume en la agricultura (WWAP2009 [2009]), será la agricultura la primera en sufrir esta reducción en forma de restricciones al riego. Ni que decir tiene que la perspectiva de poder ofrecer a los agricultores un calendario detallado con recomendaciones sobre cuándo y cómo aplicar esas restricciones tiene un valor económico inmediato. El siguiente paso es más ambicioso: si anteriormente se ofrecía al agricultor una respuesta a un problema inmediato (las restricciones), a continuación se plantea una forma de ayudarle a anticiparse a la restricción. Si fuera posible prever que a lo largo de los próximos años la disponibilidad de agua para el riego fuera a reducirse paulatinamente, se podría aprovechar el conocimiento que adquirido acerca de los efectos del riego deficitario a largo plazo en cítricos para planificar un calendario anual de riegos que permita adaptar el cultivo a la previsible nueva situación de forma que el inevitable efecto negativo sobre el cultivo y su producción sea el menor posible. En este sentido, el modelo basado en los efectos a largo plazo tomaba como punto de partida la observación de un efecto oscilante de las restricciones hídricas sobre la producción; sin embargo, en nuestro modelo este efecto oscilante se veía complementado con una pendiente creciente. Según Kozlowski y Pallardy [2002], los cultivos arbóreos que han sido expuestos con anterioridad a un estrés hídrico frecuente tienden a sufrir menores efectos en sequías posteriores. En Bréda y otros [2006] se estudia el efecto que provocaban las sequías severas sobre cultivos arbóreos silvestres (bosques) a largo plazo, y concluyó que si bien las sequías extremas podían provocar daños irreversibles y con frecuencia fatales en el cultivo, niveles de estrés hídrico moderados permitían a los árboles adaptarse a las nuevas condiciones y encontrar la estabilidad en ellas. Dickson y Tomlinson [1996] estudiaron la forma en la que las hojas, tallos y raíces se adaptan en el caso de los robles para evitar los efectos de la sequía mientras sea posible y cómo actúan para paliar dichos efectos cuando son inevitables. Es por eso que parece necesario adaptar el modelo anterior a estas propiedades de los cultivos arbóreos, buscando una estabilidad a largo plazo. La figura 7.1 muestra la comparación de las curvas de adaptación a un estrés hídrico muy leve, del 10% constante a lo largo de diez años, siguiendo dos modelos: el primero de ellos es el ya visto, en el que los efectos del estrés siguen creciendo a Modelado, simulación y control del uso del agua en la agricultura 131 Figura 7.1: Dos modelos diferentes para la pérdida de producción a largo plazo para un estrés hídrico del 10%. lo largo del tiempo con una función rampa; el segundo modelo tiene esa rampa atenuada por la exponencial negativa, lo que implica que a largo plazo el efecto creciente de la rampa queda anulado. Puede verse cómo en el modelo original, incluso con un nivel de estrés tan leve como el 10%, la pérdida de la producción -si bien se observa un comportamiento oscilante en ellasigue creciendo indefinidamente, mientras que en el segundo modelo parece estabilizarse en torno a una pérdida de producción del 6%. En este capítulo se trabajará en la búsqueda de técnicas que permitan la planificación del regadío a medio y largo plazo. En primer lugar, a partir del modelo de planificación de riego deficitario en cultivo de cítricos en función del periodo fenológico, se planteará una planificación óptima. Finalmente, se buscará aunar ambas estrategias para aprovechar nuestro conocimiento de la productividad del agua a corto plazo para mejorar la respuesta 132 Capítulo 7. Planificación de estrategias de riego del cultivo a largo plazo. Esto puede presentar ciertas complicaciones en tanto en cuanto los modelos obtenidos para la planificación en una campaña se obtuvieron a partir de ensayos de Riego Deficitario Controlado y los obtenidos para la planificación a largo plazo se obtuvieron a partir de ensayos con Riego Deficitario Sostenido. 7.2. Material y métodos 7.2.1. Planificación a largo plazo Tal y como vimos en el apartado 7.1, el modelo obtenido anteriormente no resulta útil en tanto en cuanto no se pueden planificar unas restricciones a largo plazo si la conclusión inevitable de esas restricciones es la muerte de la planta. Además, por los argumentos ya expuestos, resulta razonable pensar que el árbol, organismo vivo, puede adaptarse a las nuevas condiciones climatológicas reduciendo sus necesidades (y con ellas, su producción frutícola). El nuevo modelo propuesto, por consiguiente, sigue una ley muy similar al anterior, como puede verse en la ecuación 7.1: un modelo en variables separables xytdonde la función g1(x) es un polinomio de segundo grado y g2(t) incluye un coseno atenuado con una exponencial negativa; la diferencia es que en este caso no existe componente de rampa. L(x,t) = g1(x)g2(t) g1(x) = a1x2+a2x+a3 g2(t) = b1+e−t/τ [b3cos(ωt)+b2] (7.1) Igual que en el modelo original, hay que establecer una serie de condiciones de contorno: se asume que no se produce pérdida de producción en el momento de inicio del experimento o si el estrés es nulo. Además, se establece que la pérdida de producción en régimen permanente tiene un valor de 1 cuando la restricción es absoluta; este valor actuará como máximo de pérdida de la producción. Modelado, simulación y control del uso del agua en la agricultura 133 L(x= 0,t)=0,∀t→a3= 0 L(x,t =0)=0,∀x→b1+b2+b3= 0 →b2=−b3−b1 l´ ım t→∞L(x,t) = Lr.p.(x) = g1(x)b1→Lrp(x=1)=1→b1=1 a1+a2 (7.2) Además, como las medidas se realizan año a año, tiene sentido postular que los extremos locales (mínimos y máximos debido a la componentes oscilatoria) se encuentren en t=k,ω =π+ 2nπ,k ∈Z. Llevadas estas condiciones a 7.1, se obtiene la expresión de la ecuación 7.3. L(x,t)=(a1x2+a2x){1 a1+a2+e−t/τ [b3cos(πt)−b3−1 a1+a2]} (7.3) Considerando que la función fue creada para estimar la pérdida de producción durante un año, no tiene sentido evaluar la función L(x,t) en tiempo continuo. Llevando a cabo un muestreo con t= 1 año, se llega a la ecuación 7.4. L(x,k)=(a1x2+a2x){1 a1+a2+e−k/τ [b3(−1)k−b3−1 a1+a2]},k ∈Z(7.4) Esta función permite estimar la pérdida de fruto en un año k cuando se aplica el mismo nivel de estrés en los anteriores (k-1) años. Sin embargo, esto hace imposible predecir el resultado de variar el nivel de estrés hídrico. Para resolver este problema, es preciso transformar esta ecuación general en una ley de recurrencia, donde la pérdida en el año n sea una función de las pérdidas en años anteriores. Este paso no es inmediato. Al no tratarse de un sistema lineal e invariante en el tiempo, no tiene sentido emplear la transformada Z. La forma más simple, tal y como se hizo en el capítulo 5, es trabajar con la hipótesis de baja variación de g1(x). Por ejemplo, si asumimos que g1(xk)≈g1(xk−2), se puede poner la