基于改进遗传算法的温室环境动态优化控制

晋 春， 毛罕平， 马国鑫， 王奇瑞， 石 强

(1. 江苏大学 农业工程学院， 江苏 镇江 212013； 2. 江苏科技大学 电子信息学院， 江苏 镇江 212003)

温室通过调控室内环境以促进作物生长发育，提升作物产量和品质，从而获得高经济效益，但是，运行成本一直是制约其发展的主要因素.在保证作物高产条件下降低温室环境调控成本，以寻求经济效益最大化，一直是国内外学者长期关注的热点问题[1].本质上，这是个动态优化问题，也称之为最优控制问题，即通过建立性能目标函数，在一定的约束条件下寻求使性能目标最优的控制策略.

温室环境的最优控制，最早可以追溯到20世纪80年代，文献[2]首次建立以经济效益为最优目标的温室环境控制问题，并利用Pontryagin极大值理论(Pontryagin′s maximum principle，PMP)求解该问题.文献[3]拓展了经济效益最优目标的温室环境控制问题，并进一步分析了其PMP最优条件和求解方法，仿真结果表明比传统方法可提高经济效益多达15%.文献[4-5]利用伪谱法工具箱PROPT计算，实现玫瑰温室中输入能源最低下的最优控制方案.文献[6-7]利用工具箱PROPT实现中国日光温室内自适应闭环控制、带LED补光的双闭环最优控制，在减少能耗成本、提高作物经济效益、在线可实施性等方面都取得不错的结果.文献[8]运用差分进化(differential evolution，DE)算法进行了短时间尺度下经济效益目标最优的温室环境控制仿真研究.

动态优化问题的求解方法可以归纳为两大类[9]：间接法和直接法.间接法主要通过Pontryagin极大值理论将最优控制问题转化为两点边界问题，进行解析求解；直接法则通过离散化方法，将最优控制问题近似转化为非线性规划(nonlinear programming，NLP)问题，进行数值近似求解.常用的直接法主要有控制向量参数化法、有限元正交配置法、伪谱法、智能优化方法等.一般来说，间接法求解精度高，但是依赖于目标问题的梯度信息，最优条件计算复杂，只适用低阶的或者复杂性较低的最优控制问题；而直接法对初值的精度要求更低，其收敛域更宽广，计算效率高，是当前求解最优控制问题的主流方法[10].智能优化算法，因其不依赖于问题的梯度信息，并具有全局优化、计算简单、易编程实现等优点，在最优控制问题的求解上也引起了广泛关注[11]，同时也在温室环境最优控制上得到应用[12-13].该类方法的思想是采用分段常量近似表示控制变量，从而形成NLP问题，再利用诸如遗传算法(genetic algorithm，GA)、DE、粒子群算法(particle swarm optimization，PSO)等智能算法直接求解.其中，GA方法成熟，适用于混合整数变量等复杂优化问题.

随着温室环境智能化控制需求的不断提升，对温室环境最优控制问题求解效率和控制精度提出更高的要求.温室环境系统复杂，调控设备具有开关切换、连续与非连续控制等混合特性，导致常规的优化算法计算复杂、求解效率低.此外，其控制约束较多，存在许多相互制约作用，如CO2供给与通风、通风与加热等之间相互制约作用，致使常规的优化算法控制精度较差，容易引起设备操作不当、相互冲突等问题.因此文中运用改进遗传算法(improved genetic algorithm，IGA)，以实现温室环境动态优化问题的有效、实用求解.引入一种精确惩罚函数，以解决GA不能直接处理状态变量路径约束的问题；引入工程经验规则、精英保留、多种群进化等改进策略，以解决标准遗传算法(standard genetic algorithm，SGA)易出现未成熟收敛和后期收敛速度慢的问题，并获得满足实际控制需求的控制方案.

1 问题描述

1.1温室环境系统模型

当前，温室环境系统建模主要以机理模型的研究为主[14].考虑到文中的控制目标函数含作物的经济效益指标，这里采用含有作物生长变量的van Henten模型[3]，表达为

(1)

(2)

(3)

(4)

式中：Xd是温室作物的干物质量，kg·m-2；XT是室内空气温度，℃；Xh是室内空气绝对湿度，kg·m-3；Xc是室内CO2质量浓度，kg· m-3；Uq是加热速率，W·m-2；Uc是CO2供给率，kg· m-2·s-1；Uv是特定通风率，m·s-1；Vrad是室外太阳辐射强度，W·m-2；VT是室外空气温度，℃；Vc是室外CO2质量浓度，kg·m-3；Vh是室外空气湿度，kg·m-3.各模型参数的定义与数值, 详见文献[3]和文献[6].

在van Henten模型的基础上，定义nx维的状态变量x∈Rnx(nx=4)、nu维的控制变量u∈Rnu(nu=3)、nd维的外部天气输入量d∈Rnd(nd=5)，分别为

x=[XdXTXhXc]T，

(5)

u=[uqucuv]T，

(6)

d=[VradVTVcVhVw]T，

(7)

式中：Vw是室外风速，m·s-1；控制变量u∈Rnu由加热操作uq、通风窗操作uv和CO2供给操作uc构成.uq设置为开关控制量，uq∈[0,1]，则Uq=uq·Qheat,q，Qheat,q为恒定加热功率，150 W·m-2.uc设为开关控制量，uc∈[0,1]，则Uc=ucφinj,c，φinj,c为恒定CO2供给速率,1.2×10-6kg·m-2·s-1.

通风窗操作uv设置为连续控制变量，uv∈(0,2).它是通风窗操作的综合指令，决定了迎风窗开度uws和背风窗开度uls：

(8)

(9)

这样，特定通风率Uv再由uws和uls计算出：

(10)

那么，温室环境系统模型可简化表示为

(11)

式中：f(x,u,d,t)代表温室环境系统的系列状态方程式(1)-(4)；t是时间变量.

1.2 温室环境动态优化问题

温室环境控制的本质是为作物提供适宜的室内生长环境，以获得高经济效益，但是能耗成本也不容忽视.这是个多目标的动态优化问题.在预知作物销售价格、能耗价格系数时，以经济效益，即作物经济收入减去能耗总成本，作为单目标优化函数(性能指标J)，描述为

(12)

Φ(Xd(tf))=cpr1+cpr2Xd(tf)，

(13)

L(x,u,t)=cquq+cco2uc，

(14)

式中：t0为起始时刻；tf为终止时刻；Xd(tf)是终止时刻的作物干物质值.作物经济收入函数Φ(Xd(tf))是作物市场价格cpr1、cpr2和干物质Xd(tf)的拟合函数；能耗成本函数L(x,u,t)主要考虑加热uq的成本、CO2供给uc的成本，不计通风成本[3].

为保证室内环境适宜作物生长，避免出现极端恶劣情况，在整个时域内温室环境因子必须控制适当范围内，这在动态优化问题中被定义为状态变量x的路径约束，形如：

xi,min≤xi≤xi,max，i∈nx，

(15)

式中：xi,max、xi,min表征状态变量x约束范围的上、下界值.对于文中3个状态变量XT、XRh、Xc，对应的边界值参如表1所示.表1中XRh是室内空气的相对湿度，由空气绝对湿度Xh与饱和水气压Xh,sat的比值转换得到；Xc进行了量纲一处理.

表1 变量约束范围的边界值

根据前面的定义，控制变量u∈Rnu中uq∈[0,1]、uc∈[0,1]都为开关控制量(整数型).uv∈(0,2)为连续控制量.控制变量的约束条件统一表达为

uj,min≤uj≤uj,max，j∈nu，

(16)

式中：uj,max和uj,min表征控制量u约束范围的上、下界，也列于表1中.

进而，经济效益最优目标下的温室环境动态优化问题可以表述为问题1(P1)：在给定温室环境系统动态方程组式(1)-(4)、初始状态值x(t0)=x0时，满足控制变量约束条件式(13)和状态变量约束条件式(12)下，寻求最优控制输入量u*(t)，以使性能指标J最大化.综合表达为

xi,min≤xi(t)≤xi,max,i∈nx,

uj,min≤uj(t)≤uj,max,j∈nu.

(17)

2 算法设计与实现

遗传算法求解动态优化问题的基本思想是首先将控制变量在控制时域上离散化，转化为有限维参数的输入向量组，使动态优化问题转化为有限维参数的NLP问题，再利用模拟生物基因进化的遗传操作进行问题解的寻优.

2.1 控制变量离散化

首先对控制时域[t0,tf]分段化，形成一组有限维的时间网格节点.这样在各时间节点上，控制变量用一类基函数近似表示，被离散化成一组有限维的输入参数向量，从而将动态优化问题转化为有限维参数的NLP问题.根据基函数的不同，控制变量有多种近似策略，如分段常数、分段线性、分段样条函数等.对于遗传算法求解动态优化问题，常采用分段常数近似法，如图1所示.

图1 分段常数近似法

具体实现过程如下：

首先，将整个控制时域[t0,tf]等分成Nc个子区间[tk-1，tk],k=1,2,…,Nc，Nc记为时间网格数.时间节点tk满足：

0=t0≤t1≤…≤tNc-1≤tNc=tf，

(18)

在Nc确定后，tk为固定的时间切换点.

假设控制变量u∈Rnu的第j个分量为uj(t)(j=1,2,…,nu)，根据分段常数近似法，uj(t)可表示为

(19)

记：参数向量组σ∈Rnu×Nc，σ=[σq，σv，σc]，则控制时域[t0，tf]上连续时域变量u∈Rnu就由离散时间节点tk(k=0,1,…,Nc-1)上σ确定，可表示为u(tk|σ).因此，问题1(P1)所表示的动态优化问题，转换为有限维参数σ的非线性规划问题2(P2)：

P2: maxJ(u)=Φ(Xd(tf|σ))-L(x(t|σ),σ,t)，

xi,min≤xi(t|σ)≤xi,max,xi(t|σ)∈RNc,i∈nx,

(20)

其中：

Φ(Xd(tf|σ))=cpr1+cpr2Xd(tNc|σ) ，

(21)

(22)

式中:x(t|σ)表示在时刻t由参数向量组σ计算所得的状态变量x值.

2.2 改进遗传算法

有限维参数σ的非线性规划问题2(P2)，可采用GA等智能算法进行求解.遗传算法是一种借鉴生物界自然选择和遗传机制的随机全局搜索和优化算法，不依赖于问题的梯度信息，通过随机选择、交叉和变异操作，使种群不断繁衍进化，求得问题的最优解.编/解码、选择、交叉、变异和适应度函数是遗传算法的主要组成部分.尽管遗传算法有许多优点，但早熟收敛和后期搜索迟钝是标准遗传算法(SGA)主要存在的2个问题，通常对选择、交叉、变异3个操作进行改进，以提高算法的有效性.此外，考虑到温室环境控制的特殊性，引入工程经验规则，修正种群个体的基因(码串)，以获得满足工程要求的控制方案.

1) 编/解码.遗传算法通常使用固定长度的二进制符号串来表示种群中的个体.对于文中控制变量uq∈[0,1]、uc∈[0,1]，都为开关控制量，编码时可直接采用一位的二进制数表示.群体中的个体数量Nvars，也就是控制变量u∈Rnu离散化后的参数向量组σ∈Rnu×Nc长度，可计算为

Nvars=nuNc.

(23)

为了获得满足工程要求的解，需要对随机产生的种群个体基因(码串)进行修正，这里引入工程经验规则：CO2供给时停止通风，以及夜间不供给CO2.

2) 选择.选择是根据适应度函数选择适应性强的个体来繁殖下一代的过程.选择操作的目的就是为了避免有用遗传信息的丢失，提高全局收敛性和计算效率.SGA采用轮盘赌选择方法，令fi表示群体中第i个体的适应度值，∑fi表示群体的适应度函数值的总和，则它被选中的概率就为其适应度值所占的份额fi/∑fi.该选择方法的误差较大，有时甚至连适应度较高的个体也选不上.有效的改进方法之一是采用有条件的精英保留策略，直接将最佳的个体传递到下一代.另外，多种群迁移策略也是一种有效方法，通过多个群体并行进化、子群体间交换优良个体，增加群体的多样性，加快较好个体在群体中传播，从而提高收敛速度和解的精度.

3) 交叉.交叉是把 2个父代个体的部分结构加以交换重组而生成新个体.交叉的目的是为了能产生下一代新的个体，提高算法的全局搜索能力.交叉操作是遗传算法区别于其他进化算法的重要特征，常用单点交叉、双点交叉、均匀交叉、算法交叉等方式.交叉率pc决定交叉操作被使用的频度，pc越大，产生新个体的速度越快，但优良个体的结构被很快破坏；pc过小会使搜索过程缓慢.pc影响算法的性能和收敛性，针对不同的优化问题，需要反复试验来确定，一般取值范围为0.40～0.99.

4) 变异.变异是以较小概率对某些个体的某些位进行改变, 进而产生新的个体，是产生新的个体的辅助方法，但是必不可少的步骤，可以避免搜索过程陷于局部最优.相对于交叉运算决定了算法的全局搜索能力，变异运算决定了算法的局部搜索能力，维持群体多样性的同时，防止出现早熟现象.同样，变异率pm影响着算法的性能和收敛性，也需要多次试验来确定，一般取为10-4～10-1.

5) 适应度函数.适应度是度量群体中个体达到或者有助于找到最优解的优良程度，遗传算法在搜索过程中基本不利用外部信息，仅以个体的适应度来搜索.计算个体适应度的函数即为适应度函数，通常由目标函数变换而成的.对于Matlab遗传算法工具箱GADS，适应度函数要求非负、最小化，因此取性能指标J的倒数为适应度函数，即f=1/J.

2.3 状态变量路径约束处理

对于形如式(15)的状态变量不等式约束，GA不能直接处理.一般利用静态罚函数方法来间接处理，即在原目标函数基础上增加惩罚项，构成新的目标函数.这里，采用精准罚函数[15]：

(24)

式中：Gi为状态变量的第i个不等式约束条件；ρ为惩罚系数.对于双边界约束的式(15)，需要转变为2个Gi≤0式.

对应于文中3个状态变量XT、Xh和Xc，在时间节点tk上的惩罚值分别记为P(XT(tk|σ))、P(Xh(tk|σ))和P(Xc(tk|σ)).增加罚函数项后，式(19)中能耗成本函数L(x(t|σ),σ,t)更新为

(25)

2.4 算法实现

考虑到Matlab自带遗传算法工具箱GADS的通用、便捷、易编程等特点，用户只需要设置各类参数就能实现遗传算法的优化计算功能.文中将利用工具箱GADS完成算法的实现，精英保留策略、多种群迁移等改进策略也通过相关的参数设置来实现.运用遗传算法实现温室环境最优问题的求解时，主程序流程大致由控制变量离散化、GA参数设置、GA命令调用、结果输出与显示等部分组成.以改进的遗传算法(IGA)为例，详细实现步骤如下：

step1控制变量离散化.设定控制时域[t0，tf]、时域区间分段数Nc.确定控制变量u∈Rnu，采用分段常数近似法，产生GA输入量σ∈Rnu×Nc.按照σq、σc、σv顺序组合成一维数组形式.根据表1中各控制变量u∈Rnu约束范围，设定输入量σ的取值范围[LB,UB].

控制变量uq∈[0,1]、uc∈[0,1]，都为开关控制量，需要设定对应的参数组σq、σc为整数型.因此，指定GA控制参数“IntCon”为[1，2Nc]，即可将输入量σ中的前1到2Nc个分量设定为整数取值.

step2设定GA控制参数“options”，主要包含群体大小“PopulationSize”、迭代次数“Generations”、收敛精度“TolFun”、交叉率“CrossoverFraction”等；其他参数都采用默认值.

通过参数“EliteCount”设置为一个正整数，可以实现精英保留策略的改进方法.默认值为2，表示每代中选择2个优良个体传递到下一代.

设置群体大小“PopulationSize”为一个长度大于1的向量，则可实现多种群并行进化策略，子种群的数量等于向量的长度，每个子种群的大小是向量的对应项值.子种群间个体迁移可以通过设置方向“Direction”、间隔“Interval”、百分比“Fraction”来进行控制.

step3调用工具箱GADS的GA命令函数，完成优化计算.GA命令调用的基本格式为

[x,fval]=ga(@objfun,Nvars,[], [],[],[],

LB,UB, IntCon, options)，

其中：x,fval为最优解及其对应的适应度值；objfun为种群个体的适应度计算函数.对于工具箱GADS，只需要执行上述GA命令就可实现遗传算法的全部操作.而objfun函数是用户程序设计的核心，其具体实现步骤为① 读取遗传操作后的种群，运用工程经验规则，修正种群编码；进行解码、生成参数向量组σ.② 读取室外气候数据d，通过插值方法配置在时间节点tk(k= 0,1,…,Nc)上.③ 基于温室环境系统模型，求解状态变量x∈Rnx×Nc.设定t0时刻的状态变量初始值x0，根据微分方程式(1)-(4)，采用四/五阶变步长龙格-库塔法，顺序求解出各时间节点tk(k= 1,2,…,Nc)上状态变量值x.④ 适应度计算. 基于式(18)、(22)、(12)，分别计算作物经济收入Φ、能耗成本L(含惩罚项)和经济效益目标值J.将经济效益最大化的目标值J进行倒数变换，返回适应度值f.

step4输出最优解，计算并显示相关结果.

3 优化仿真试验

3.1 仿真试验条件

计算机配置为酷睿双核i5-4460，主频3.2 GHz，内存12 G.仿真试验平台采用Matlab2016a软件，GA算法实现采用Matlab自带的遗传算法工具箱GADS.

用于算法仿真的室外天气数据，来自室外气象站(江苏省镇江市江苏大学农业工程学院温室)记录的实际数据.控制时域取3 d，t0=0 h，tf=72 h.根据温室系统模型需求，选取较冷的3 d室外天气数值(2018年1月29日到31日)，采样时间为1 min，如图2所示.室外CO2质量浓度取恒定值360×10-6.

图2 3 d室外天气数值

3.2 仿真结果与分析

针对给定的温室环境动态优化问题2(P2)的求解，在SGA基础上，逐步增加处理状态变量路径约束的惩罚项、带有工程约束规则和多种群进化的改进策略，实现未含有惩罚项的标准遗传算法SGA(no)、含有惩罚项的标准遗传算法SGA(penalty)和含有惩罚项的改进遗传算法IGA(penalty)3种算法的对比.在设定的控制时间3 d下，时间网格数Nc都取为144个，即每个控制变量u∈Rnu离散为144个参数的向量组，则计算出GA输入量的个数Nvars为432个.

GA控制参数设置影响到算法性能.为了获得良好的计算效果，并提高程序的计算速率，文中通过对不同参数组合重复试验，确定种群大小取100，迭代次数取200次，收敛精度取1×10-4.对于采用多种群并行进化的改进遗传算法，设置2个种群，大小为50，精英保留个数为10个，迁移率为0.1，其他参数采用GA工具箱默认值.

在图2所示的室外天气作用下，状态变量的求解结果如图3所示.

图3 3种方法的状态变量优化轨迹

图3b-d中水平线虚线表征表1中状态变量约束的边界值.可见，SGA(no) 优化的状态量轨迹明显出现超出规定边界值的现象，不满足作物生长对环境值的严格要求，说明GA本身对状态变量路径约束条件的处理能力非常弱.而增加惩罚项后，SGA(penalty)和IGA(penalty) 的优化结果始终保持在规定的界限范围内，表明惩罚函数大大提高了GA处理路径约束条件的能力，提升了GA算法的有效性、实用性.从图3b、d中室内温度XT、室内CO2质量浓度Xc的优化轨迹来看，SGA(penalty) 的优化结果没有明显的规律性，而IGA(penalty) 的优化结果，特别是Xc，在2个白天时段(9～18 h、57～64 h)维持较高的数值，与室外太阳辐射强度Vrad关联性很大，符合温室环境控制的实际物理规律.通常，室外光照强时，需要较高的室内温度、CO2质量浓度，这将有助于作物光合作用，从而促进作物生长发育.反之光照弱时，作物光合作用也较弱，作物生长缓慢，无须维持较高的室内温度、CO2质量浓度.

图3a表明了优化的结果也符合作物生长发育规律.在选定的3 d中，9～16 h、57～64 h 2个白天时段内较高的室外太阳辐射强度、室内温度、室内CO2质量浓度将促使作物生长，因此干物质Xd增长较快.而在32～41 h白天时段，室外光强低，由于所优化出室内温度、室内CO2质量浓度也较低，致使光合作用产物较低，Xd增长非常小.图3c主要显示出空气相对湿度的变化状况，通常空气温度较高时，相对湿度就较低.

各方法的控制变量优化结果如图4所示.

图4 3种方法的控制变量优化轨迹

总体来看，作为随机搜索算法的GA，优化结果的波动性较大.由于性能指标中不包含通风成本，致使图4b中各方法的通风操作随机性都较大.图4a中IGA(penalty) 的加热uq动作较多，是为了维持较高的室内温度所致.图4c中，相比于其他2种方法，IGA(penalty) 的CO2供给uc动作最少，只需要维持2个白天时段(9～18 h、57～64 h)内CO2质量浓度较高即可，另外与图4b中uv在时序上不冲突.而其他2种方法优化的uv和uc冲突情况很频繁，严重削弱了CO2供给的作用效果，在经济上是不可取的，也不符合实际温室控制中的物理规律.可见，带有工程约束规则的改进遗传算法，有效提高了控制品质，可获得满足工程要求的控制方案.

进一步优化结果中几个关键的性能指标如表2所示.为避免算法单次运行的差异性，统计了算法运行10次的结果.表2中迭代次数是算法满足终止条件时的迭代运算次数，表征算法的收敛速度. 在SGA(no)多次仿真结果中，状态变量轨迹多次超出设定界限，温度值甚至出现低于0 ℃，导致模型计算出的总干物质Xd为负值，失去比较意义.因此，这里仅比较SGA(penalty)和IGA(penalty)的优化结果.可以看出，相对于SGA(penalty)，IGA(penalty)优势较为明显，迭代次数降低接近20次，体现出多种群并行进化等改进策略的作用效果，有效提升了算法的收敛速度.此外， 虽然IGA(penalty) 的总干物质Xd平均值下降了2.51%，但是能耗成本L、经济利润J全部指标值占有优势，能耗成本L平均值降低了约12.62%，最大可降低41.36%；经济利润J平均值略提高了0.60%，最大能够提高1.34%.可见，带有改进策略的IGA(penalty)算法有效提高了性能指标值.

表2 2种方法的性能指标对比

综合来看，针对经济效益最优目标的温室环境控制问题，采用惩罚函数的遗传算法求解方法可以有效满足状态变量路径约束条件，但是控制变量的解还会违背了温室控制实际的物理规律，说明GA算法本身在温室环境动态优化控制问题的求解上，有效性、实用性不能满足工程应用要求，这可能也是智能算法应用于温室环境优化控制上相关研究文献较少的一个原因.而通过引入工程经验规则、精英保留、多种群并行进化等改进措施，IGA(penalty)算法优化效果提升明显，算法性能和控制品质得以有效提升，控制变量的优化结果更符合实际的温室控制物理规律，控制方案可行性更高.

4 结 论

针对一种整数变量控制下经济效益最优目标的温室环境控制问题，在标准遗传算法的层面上，采用惩罚函数以处理状态变量路径约束条件；引入整数变量设置和工程经验规则来修正遗传编码，以构建满足实际控制需求的解空间；采用精英保留、多种群并行进化等改进措施以提升算法性能.仿真试验表明，相比标准遗传算法，所提方法能够有效提升收敛速度、提高控制品质，能耗成本最大可降低41.36%，经济利润最大能提高1.34%，显示出其潜在的应用价值.