CN104899668A - 基于污染税的流域污染负荷分配二层决策优化方法 - Google Patents

Abstract

本发明提出了基于污染税的流域污染负荷分配二层决策优化方法，首先，系统获取水流量，并将其转换为水流量的确定性变量；其次，系统获取河流中溶解氧浓度关系及生化需氧量浓度关系；再其次，系统根据水流量的确定性变量、溶解氧浓度关系及生化需氧量浓度关系构建下层模型；然后，系统根据污染税收标准制定上层模型；再然后，系统根据下层模型及上层模型构建基于污染税的流域污染负荷分配二层决策模型；最后，系统对基于污染税的流域污染负荷分配二层决策模型进行求解，得到最优解。适用于基于污染税的流域污染负荷分配问题。

Description

基于污染税的流域污染负荷分配二层决策优化方法

技术领域

本发明涉及流域水质目标管理，特别涉及基于污染税的流域污染负荷分配决策方法。

背景技术

二层决策优化方法是一类针对具有二层递阶结构问题的优化方法。其主要研究具有两个层次系统的规划与管理问题。分为上层和下层，上下两层问题都具有各自的决策变量、约束条件和目标函数。这种方法区别于其他优化方法的一个主要方面是不再是一个决策者进行规划与决策，而是有上下两层决策者，上层决策者只是通过自己的决策去指导下层决策者，并不直接干涉下层的决策；而下层的决策将上层的决策作为一个影响条件，其在自己的可行域进行自由决策，这就导致上层决策者在做任何决策的时候必须考虑到下层决策者可能做出的反馈，以规避下层决策者的决策可能带来的不利影响。

在水污染严重的今天，合理的流域水质管理显得至关重要。目前，为进一步治理各流域水环境，我国环保部制定了一系列的环境经济政策，包括环境税费、排污权交易等。

其中，污染税费为主要方式。其旨在国家以立法的形式明确规定污染者要按其排放的污染物种类、浓度和数量缴纳一定的费用。通过征收污染税，使经济活动所造成的环境损失由其本身承担，使外部经济内在化；此外，税金的收入可用作改善环境项目的拨款。污染税费制定在限制各流域排放者的排污行为的同时，也一定程度的增加了地方环保部门的水污染治理工作难度。具体表现为：①如何设计一套科学的收费体系，可以使某一流域和地区各排污者实现特定水质标准所需的成本相同；②如何制定一套合理的流域污染负荷控制技术，能够指导地方环保部门在污染税费政策下进行水污染防治工作。因此，本发明考虑了基于污染税的流域污染负荷分配问题，旨在提出一种新的水环境污染分配方案。

在过去的几十年中，决策优化技术在水资源管理方面起到了至关重要的作用，主要相关的技术包括多目标决策优化，多属性决策优化以及多阶段决策优化等，这些决策的一个主要特点是由一个主体进行决策，没有考虑其他利益相关者的意见及反馈。

发明内容

本发明所要解决的技术问题，就是提供一种基于污染税的流域污染负荷分配二层决策优化方法以解决基于污染税的流域污染负荷分配问题。

本发明解决所述技术问题，采用的技术方案是，基于污染税的流域污染负荷分配二层决策优化方法，包括以下步骤：

步骤1、系统获取水流量，并将其转换为水流量的确定性变量；

步骤2、系统获取河流中溶解氧浓度关系及生化需氧量浓度关系；

步骤3、系统根据水流量的确定性变量、溶解氧浓度关系及生化需氧量浓度关系构建下层模型；

步骤4、系统根据污染税收标准制定上层模型；

步骤5、系统根据下层模型及上层模型构建基于污染税的流域污染负荷分配二层决策模型；

步骤6、系统对基于污染税的流域污染负荷分配二层决策模型进行求解，得到最优解。

具体的，所述步骤1中，首先，系统采用模糊随机变量对未来时间水流量进行量化，系统将水流量分为三个不同的级别:低水位(L)，中水位(M)和高水位(H)，服从概率分布(a,b,c)，三种水位的概率分别可以表示为p_L，p_M和p_H，综上，水流量为:

其次，将水流量转换为模糊随机变量，并计算其模糊期望值，具体包括：

①假设水流量可以表述成为任意一个模糊随机变量χ，含有以下数值：

${\tilde{Q}}_L=(a_L,b_L,c_L),{\tilde{Q}}_M=(a_M,b_M,c_M),{\tilde{Q}}_H=(a_H,b_H,c_H);$

其中，表示三个水流量水平，因此其概率分布如下：

$P\left(\right\{\mathrm{\ω}\∈\mathrm{\Ω}\left|\mathrm{\χ}\right(\mathrm{\ω})={\tilde{Q}}_L\})=p_L,P\left(\right\{\mathrm{\ω}\∈\mathrm{\Ω}\left|\mathrm{\χ}\right(\mathrm{\ω})={\tilde{Q}}_M\})=p_M,P\left(\right\{\mathrm{\ω}\∈\mathrm{\Ω}\left|\mathrm{\χ}\right(\mathrm{\ω})={\tilde{Q}}_H\})=p_H;$

其中，Ω表示三种水流水平及其模糊逻辑；

假设：θ(χ)是模糊随机变量χ的模糊参数，如果θ(χ)＝E^f[X|P],存在，则有：

$inf{\left(\mathrm{\θ}\right(\mathrm{\χ}\left)\right)}_{\mathrm{\α}}=E^f\left({\mathrm{inf\χ}}_{\mathrm{\α}}\right|P)=p_{1}inf{\left({\tilde{q}}_1\right)}_{\mathrm{\α}}+p_{2}inf{\left({\tilde{q}}_2\right)}_{\mathrm{\α}}+p_{3}inf{\left({\tilde{q}}_3\right)}_{\mathrm{\α}}$

${\sup \left(\mathrm{\θ}\right(\mathrm{\χ}\left)\right)}_{\mathrm{\α}}=E^f\left({\mathrm{sup\χ}}_{\mathrm{\α}}\right|P)=p_1\sup {\left({\tilde{q}}_1\right)}_{\mathrm{\α}}+p_2\sup {\left({\tilde{q}}_2\right)}_{\mathrm{\α}}+p_3\sup {\left({\tilde{q}}_3\right)}_{\mathrm{\α}};$

那么，模糊随机变量χ的模糊期望值可以表示为则有:

$E^I\[\widetilde{\stackrel{\‾}{Q}}\]=(a_E,b_E,c_E);$

其中，a_E＝inf(θ(χ))₀,b_E＝inf(θ(χ))₁＝sup(θ(χ))₁,c_E＝sup(θ(χ))₀，表示模糊随机变量χ的模糊期望值；

②模糊随机水流量被转化成为三角模糊数，系统采用的带乐观-悲观指数的期望值算子将其转化成为确定性值；

首先，定义一个乐观-悲观指数(0≤ε≤1)用于量化决策者的决策偏好，然后使用期望值算子计算该三角模糊数的期望值，并将其作为水流量的确定性变量，则有：

$E^{\mathrm{\ϵ}}\[E^I\[{\widetilde{\stackrel{\‾}{Q}}}_u\]\]=\frac{1-\mathrm{\ϵ}}{2}{a}_E+\frac{1}{2}{b}_E+\frac{\mathrm{\ϵ}}{2}{c}_E;$

其中，表示水流量的确定性变量。

具体的，所述步骤2中，系统根据Streeter-Phelps模型的结论，河流中生化需氧量浓度的随时间变化量与其自身浓度成正比，即：

$\frac{{\mathrm{dC}}_b}{dt}=-k_d{C}_b;$

其中，C_b表示生化需氧量浓度，k_d表示脱氧系数；

因此，在一个河段t，用k_d(t)表示脱氧系数，d(t)表示该河段长度，v(t)表示该河段评价水流流速，C_b(t)表示流入该河段水流的生化需氧量浓度，则流出该河段水流的生化需氧量浓度C_b(t+1)为定义该生化需氧量浓度等价于：

C_b(t+1)＝α(t)C_b(t)；

系统根据Streeter-Phelps模型的结论，河流中的溶解氧从空气中补充，其补充浓度的随时间变化量与其氧亏浓度成正比，即：

$\frac{{\mathrm{dD}}_o}{dt}=-k_d{C}_b+k_a{D}_o;$

其中，C_b表示生化需氧量浓度，D_o表示氧亏浓度，k_d表示脱氧系数，k_a表示充氧系数。与生化需氧量情况相似，在一个河段t，用k_a(t)表示充氧系数，C_b(t)，C_o(t)分别表示流入该河段水流的生化需氧量浓度，C_do表示该河段水流的饱和溶解氧浓度，则流出该河段水流的溶解氧浓度C_o(t+1)为：

$C_o(t+1)=e^{-k_a\left(t\right)\frac{d\left(t\right)}{v\left(t\right)}}{C}_o\left(t\right)-\frac{k_d\left(t\right)}{k_a\left(t\right)-k_d\left(t\right)}\[e^{-k_d\left(t\right)\frac{d\left(t\right)}{v\left(t\right)}}-e^{-k_a\left(t\right)\frac{d\left(t\right)}{v\left(t\right)}}\]C_b\left(t\right)-C_{do}\[1-e^{-k_a\left(t\right)\frac{d\left(t\right)}{v\left(t\right)}}\];$

定义：和δ(t)＝C_do[1-β(t)]，该溶解氧浓度等价于：

C_o(t+1)＝β(t)C_o(t)-γ(t)C_b(t)-δ(t)；

其中，α(t)、β(t)分别可以表示河段t生化需氧量浓度(BOD)、溶氧浓度(DO)的转换系数，γ(t)、δ(t)分别表示影响河段t溶氧浓度的生化需氧量浓度、饱溶氧浓度转换系数。

具体的，所述步骤3中，系统构下层模型包括：

(1)河流系统状态方程组构建，具体包括：

①流量状态转移方程：基于质量守恒定律以及流量连续性条件，用Q_r(t-1)，Q_r(t)分别表示在监测点t-1和t的流量，用Q_w(t)，Q_d(t)分别表示排污控制单元t的取水量和排污量，则可构建流量状态转移方程如下：

Q_r(t)＝Q_r(t-1)-Q_w(t)+Q_d(t)；

②BOD浓度状态转移方程：基于质量守恒定律、Streeter-Phelps水质模拟模型以及流量状态转移方程，用C_b(t-1)，C_b(t)表示监测点t-1和t的BOD浓度，表示排污控制单元t未处理的污水BOD浓度，以及η(t)为排污控制单元t进行削减污水BOD浓度的比例，其为下层排污者t的决策变量，则可构建BOD浓度状态转移方程如下：

$C_b\left(t\right)Q_r\left(t\right)=\mathrm{\α}(t-1)C_b(t-1)\[Q_r(t-1)-Q_w\left(t\right)\]+C_b^i\left(t\right)\[1-\mathrm{\η}\left(t\right)\]Q_d\left(t\right);$

定义变量则该状态转移方程等价于：

$C_b\left(t\right)=\[1-r\left(t\right)\]\mathrm{\α}(t-1)C_b(t-1)+r\left(t\right)C_b^i\left(t\right)\[1-\mathrm{\η}\left(t\right)\];$

③DO浓度状态转移方程：相似于BOD情形，基于质量守恒定律、Streeter-Phelps水质模拟模型以及流量状态转移方程，用C_o(t-1)，C_o(t)表示监测点t-1和t的DO浓度，表示排污控制单元t未处理的污水DO浓度，则可构建DO浓度状态转移方程如下：

$C_o\left(t\right)Q_r\left(t\right)=\[\mathrm{\β}(t-1)C_o(t-1)-\mathrm{\γ}(t-1)C_b(t-1)-\mathrm{\δ}(t-1)\]\[Q_r(t-1)-Q_w\left(t\right)\]+C_o^i\left(t\right)Q_d\left(t\right);$

同样的定义变量则该状态转移方程等价于：

$C_o\left(t\right)=\frac{Q_r(t-1)-Q_w\left(t\right)}{Q_r\left(t\right)}\[\mathrm{\β}(t-1)C_o(t-1)-\mathrm{\γ}(t-1)C_b(t-1)-\mathrm{\δ}(t-1)\]+\frac{Q_d\left(t\right)}{Q_r\left(t\right)}{C}_o^i\left(t\right);$

另外，考虑河流中饱和溶解氧浓度C_do的限制，以上DO浓度状态转移方程可改写为：

$C_o\left(t\right)=min\{\[1-r\left(t\right)\]\[\mathrm{\β}(t-1)C_o(t-1)-\mathrm{\γ}(t-1)C_b(t-1)-\mathrm{\δ}(t-1)\]+r\left(t\right)C_o^i\left(t\right),C_{do}\};$

(2)构建初始模型河流状态条件：

对于位于整个河流系统中上游第一个排污控制单元，用分别表示上游水流BOD和DO浓度，则该控制单元的模型考虑初始河流状态条件如下：

$\begin{array}{ccc}{Q}_r\left(0\right)=E^{\mathrm{\ϵ}}\[E^I\[{\widetilde{\stackrel{\‾}{Q}}}_u\]\],& C_b\left(0\right)=C_b^u,& C_o\left(0\right)=C_o^u\end{array};$

(3)构建下层模型约束条件：

①允许的BOD和DO浓度标准：

排污控制单元t排放污水到河流中，不能使得在相应检测点河流的BOD浓度(C_b(t))超过环保部门所允许的最大BOD浓度标准(B_u(t))，即：

0≤C_b(t)≤B_u(t)；

同样，不能使相应检测点河流的DO浓度(C_o(t))低于环保部门所允许的最小DO浓度标准(D_l(t))，即：

D_l(t)≤C_o(t)；

②BOD削减能力约束：

每个排污控制单元具有污水BOD浓度削减能力的限制，通常采用二级污水处理工艺。对于排污控制单元t，第一级处理后的污水BOD浓度削减率为第二级处理后的污水BOD浓度削减率为在实际的污水处理运行中，环保部门要求所以排污者的污水必须进行第一级处理；因此，排污控制单元t进行污水BOD浓度削减决策，其BOD削减率(η(t))不能够超过同时不能低于即：

${\mathrm{\η}}_t^{\min }\≤\mathrm{\η}\left(t\right)\≤{\mathrm{\η}}_t^{\max };$

(4)构建下层模型目标函数：

下层决策者的目标为最小的单位污染成本目标函数，污染成本由污水处理成本和污染税费组成；

①污水处理成本：

基于BOD削减和污水处理成本的经验曲线，本发明对于排污控制单元的污水处理费用采用幂函数形式，用k_c1(t)和k_c2(t)表示排污控制单元t的成本参数1和成本参数2，k_ce(t)和k_cb(t)表示排污控制单元t的污水和BOD削减参数，则污水处理成本C(t)可以表示为以下方程：

$C\left(t\right)=k_{c1}\left(t\right)Q_d{\left(t\right)}^{k_{ce}\left(t\right)}+k_{c2}\left(t\right)Q_d{\left(t\right)}^{k_{ce}\left(t\right)}\mathrm{\η}{\left(t\right)}^{k_{cb}\left(t\right)};$

其次，参数k_c1(t)，k_c2(t)，k_ce(t)和k_cb(t)数值以污水处理构筑物建设和运营所得数据，采用最小二乘法进行曲线拟合所得；

②污染税费：

在每个监测点，地方环保部门制定相应的污水税收标准(Al_b(t)，上层模型的决策变量)和税率(T_b)，如果检测点t所对应的河流BOD浓度(C_b(t))低于该税收标准(C_b(t)≤Al_b(t))，则排污控制单元t不需要征收污染税(P(t))；如果检测点t所对应的河流BOD浓度(C_b(t))超过该税收标准(C_b(t)>Al_b(t))，则排污者t需征收污染税P(t)＝Q_r(t)T_b[C_b(t)-Al_b(t)]。因此，排放者t的污染税费模型定义为：

$P\left(t\right)=Q_r\left(t\right)T_b\left\{\frac{|C_b\left(t\right)-{\mathrm{Al}}_b\left(t\right)|+\[C_b\left(t\right)-{\mathrm{Al}}_b\left(t\right)\]}{2}\right\};$

综上所述，排污控制单元t的目标函数，即最小化的单位污染成本如下：

$\underset{\mathrm{\η}\left(t\right)}{min}{F}_t=\frac{C\left(t\right)}{Q_d\left(t\right)}+\frac{P\left(t\right)}{Q_d\left(t\right)}.$

具体的，所述步骤4中，系统构上层模型包括：

(1)上层模型约束条件构建：

对于上层模型的决策者-地方环保部门，在进行相应污染税收标准的制定过程中，需满足允许的BOD浓度标准约束条件，环保部门所制定的污染税收标准(Al_b(t))不能超过允许的最大BOD浓度标准(B_u(t))，即：

0≤Al_b(t)≤B_u(t)；

(2)上层模型目标函数构建：

上层决策者需保证公平的污染负荷分配，因此，本发明采用各排污控制单元的单位污染成本的方差来量化水污染负荷分配的公平性，用表示各排污控制单元平均的单位污染成本，则其最小的方差定义为：

$\underset{{\mathrm{Al}}_b\left(t\right)}{min}dF=\underset{t=1}{\overset{T}{\Σ}}{(F_t-aF)}^2;$

具体的，所述步骤5中，系统根据下层模型及上层模型构建的基于污染税的流域污染负荷分配二层决策模型为：

$\underset{{\mathrm{Al}}_b\left(t\right)}{min}dF={\mathrm{\Σ}}_{t=1}^T{(F_t-aF)}^2$

具体的，所述步骤6中，假设决策排污者为T个，则在特定上层决策变量下地下层模型为T个一元的非线性规划模型；首先，系统利用极值算法计算出单一的非线性规划模型的极值，再通过下层约束条件计算出该单一的非线性规划模型的最优解，进而获得下层T个一元的非线性规划模型的全部最优解集；利用粒子群算法对上述最优解集进行计算，从而获得基于污染税的流域污染负荷分配二层决策模型的最优解。

具体的，粒子群算法的粒子更新方程如下：

P_l(τ+1)＝V_l(τ+1)+P_l(τ)；

其中，表示当前个体最优粒子，即对o＝1,2,...,τ；表示当前粒子群中最优粒子，即对l＝1,2,...,L；θ(τ)表示惯性权重，计算过程为其中θ(1)和θ(I)表示第1和第I代的惯性权重；c_ρ和c_l分别表示个体和全局最优学习因子；r_ρ和r_l表示0到1之间的随机数。

具体的，系统采用Boltzmann选择算子代替传统全局最优粒子的更新机制。

本发明的有益效果是：本发明针对河流流量这一不确定参数，采用模糊随机数进行描述，并采用两步法进行处理，即首先将其转化成为模糊数，然后采用模糊期望值算子将模糊数转化为确定值。针对二层决策模型的求解，根据在污染税政策下水污染负荷分配的二层决策结构特征以及河流系统中污染物的传递特性，提出了混合嵌套基于多智能体的动态极值算法的改进粒子群算法(PSO)的方法进行求解，从而得到合理的流域污染税费方案和污染物削减方案。

附图说明

图1为本发明基于污染税的流域污染负荷分配二层决策优化方法实施例中基于污染税的流域污染分配问题的二层决策模型结构；

图2为本发明基于污染税的流域污染负荷分配二层决策优化方法实施例中BOD和DO转化机制的概念图；

图3为本发明基于污染税的流域污染负荷分配二层决策优化方法实施例中河流系统转态转态的概念图；

图4为本发明基于污染税的流域污染负荷分配二层决策优化方法实施例中排污控制单元污水处理分析；

图5为本发明基于污染税的流域污染负荷分配二层决策优化方法实施例中针对基于污染税的流域污染负荷分配问题的二层决策建模技术流程；

图6为本发明基于污染税的流域污染负荷分配二层决策优化方法实施例中粒子结构表示；

图7为本发明基于污染税的流域污染负荷分配二层决策优化方法实施例中粒子性能的混合嵌套计算流程图；

图8为本发明基于污染税的流域污染负荷分配二层决策优化方法实施例中更新全局最优粒子的Boltzmann选择算子流程图；

图9为本发明基于污染税的流域污染负荷分配二层决策优化方法实施例中河流系统图。

具体实施方式

下面结合附图及实施例详细描述本发明的技术方案：

本发明针对基于污染税的流域污染负荷分配问题，提出了基于污染税的流域污染负荷分配二层决策优化方法，首先，系统获取水流量，并将其转换为水流量的确定性变量；其次，系统获取河流中溶解氧浓度关系及生化需氧量浓度关系；再其次，系统根据水流量的确定性变量、溶解氧浓度关系及生化需氧量浓度关系构建下层模型；然后，系统根据污染税收标准制定上层模型；再然后，系统根据下层模型及上层模型构建基于污染税的流域污染负荷分配二层决策模型；最后，系统对基于污染税的流域污染负荷分配二层决策模型进行求解，得到最优解。针对河流流量这一不确定参数，采用模糊随机数进行描述，并采用两步法进行处理，即首先将其转化成为模糊数，然后采用模糊期望值算子将模糊数转化为确定值。针对二层决策模型的求解，根据在污染税政策下水污染负荷分配的二层决策结构特征以及河流系统中污染物的传递特性，提出了混合嵌套基于多智能体的动态极值算法的改进粒子群算法(PSO)的方法进行求解。从而得到合理的流域污染税费方案和污染物削减方案。

实施例

本发明主要研究单一河流系统中污染税政策下的水污染负荷分配设计。单一河流系统由地方环保部门(EPA)以及一系列坐落于同一流域的排污控制单元(Dischargers)组成。环保部门为了衡量各排污者污染水排放情况，在其下游安置水质监测点以监测河流水质情况，并根据流域污染物对人类和生态系统健康的影响，以及流域环境生态功能分区情况，在各监测点制定相应的水质标准。在满足水质标准条件下，每个排污者被允许将部分污染物削减后的污水排放入河流中。在污染税政策下，环保部门为化解上下游排污者之间水污染负荷分配的不公平以及进一步限制各排污者的污水排放，在各个监测点制定相关污染物污染税费标准，如果某监测点监测相应污染物浓度大于该标准，环保部门将向相应的排污控制单元征收相应的污染税收。每个排污控制单元需要平衡相应的污水处理成本和污染税费比例，决策污水处理的削减比例以实现最小成本达到相应水质标准。

一般情况下，即使地方环保局和各排污控制单元是独立的决策主体，地方环保局作为领导者具有优先决策权，所以各排污控制单元是从属者。因此，基于污染税的流域污染负荷分配问题可以抽象成为一个合作型二层多从属者规划问题，整个过程可以分为两个环节：污染税制定环节和污染物削减环节：在污染税制度环节，地方环保部门(领导者)制定T个监测点的污染税收标准以获得最小的所有排污控制单元的单位排污成本的偏差。在污染物削减环节：各排污控制单元(从属者)依据环保部门制定的税收标准以及上游水质情况依次决策自身污水削减比例以获得最小的单位排污成本，如图1所示。

为解决上述基于污染税的流域污染负荷分配问题，提出如下合作型二层多从属者决策建模技术，具体步骤如下：

1.不确定变量的量化与处理。

在水质量管理问题中需要考虑不确定问题已经得到广泛的认知，因为不确定存在于各种不同的系统组成部分和各种预决策环境当中。因此，现实水质量管理决策中的内部复杂以及随机性使其不能单单的用传统的确定性优化方法来考虑与解决。因此，如要进一步对决策中的不确定变量进行量化和处理。

(1)量化过程：在污染负荷分配决策中，对上游水流量的量化十分重要，然而，上游水量是一个典型的不确定性，其大小受到多种因素影响，如温度、降雨量、蒸发量等。因此，本发明采用模糊随机变量对未来时间水流量进行量化，即可以采取相应的水资源配置决策由于河流水流量的变化，未来的水供应是不确定的，水流量一般可以分为三个不同的级别:低水位(L)，中水位(M)和高水位(H)，服从概率分布(a,b,c)，三种水位的概率分别可以表示为p_L，p_M和p_H。然而，对于每一个水位，区域水资源管理局很难给定一个确切的值，因此需要用模糊概念来描述它。为了搜集每一个水位的数据，针对水位k，找相应的水文专家进行数据采集，他们通常用口头表述“水流量在a_k和c_k之间，并且有一个最有可能的值b_k。”其中，k表示水位(t＝L表示低水位，t＝M表示中水位，t＝H表示高水位)。值得一提的是，这些专家给定的口头描述是通过多次的观察的数据获得。因此在每一个水位上的水流量可以用三角模糊数描述，综上，水流量(可供的分配水量)可以考虑成为一个模糊随机变量。

(2)处理过程：一般地，由于模糊随机变量的存在我们很难获得模型的解。本发明提出一种两步转化方法，将模糊随机变量转化成为确定型变量。第一步：将模糊随机变量转化成为模糊数；第二步：运用带乐观-悲观指数的期望值算子将模糊数转化成为确定型变量。关于上游水流量具体的转换步骤如下：

①假设水流量可以表述成为任意一个模糊随机变量χ，含有以下数值：

${\tilde{Q}}_L=(a_L,b_L,c_L),{\tilde{Q}}_M=(a_M,b_M,c_M),{\tilde{Q}}_H=(a_H,b_H,c_H);$

其中，表示三个水流量水平，因此其概率分布如下：

$P\left(\right\{\mathrm{\ω}\∈\mathrm{\Ω}\left|\mathrm{\χ}\right(\mathrm{\ω})={\tilde{Q}}_L\})=p_L,P\left(\right\{\mathrm{\ω}\∈\mathrm{\Ω}\left|\mathrm{\χ}\right(\mathrm{\ω})={\tilde{Q}}_M\})=p_M,P\left(\right\{\mathrm{\ω}\∈\mathrm{\Ω}\left|\mathrm{\χ}\right(\mathrm{\ω})={\tilde{Q}}_H\})=p_H;$

其中，Ω表示三种水流水平及其模糊逻辑；

假设：θ(χ)是模糊随机变量χ的模糊参数，如果θ(χ)＝E^f[X|P],存在，则有：

$inf{\left(\mathrm{\θ}\right(\mathrm{\χ}\left)\right)}_{\mathrm{\α}}=E^f\left({\mathrm{inf\χ}}_{\mathrm{\α}}\right|P)=p_{1}inf{\left({\tilde{q}}_1\right)}_{\mathrm{\α}}+p_{2}inf{\left({\tilde{q}}_2\right)}_{\mathrm{\α}}+p_{3}inf{\left({\tilde{q}}_3\right)}_{\mathrm{\α}}$

$\sup {\left(\mathrm{\θ}\right(\mathrm{\χ}\left)\right)}_{\mathrm{\α}}=E^f\left({\mathrm{sup\χ}}_{\mathrm{\α}}\right|P)=p_1\sup {\left({\tilde{q}}_1\right)}_{\mathrm{\α}}+p_2\sup {\left({\tilde{q}}_2\right)}_{\mathrm{\α}}+p_3\sup {\left({\tilde{q}}_3\right)}_{\mathrm{\α}};$

那么，模糊随机变量χ的模糊期望值可以表示为则有:

$E^I\[\widetilde{\stackrel{\‾}{Q}}\]=(a_E,b_E,c_E);$

其中，a_E＝inf(θ(χ))₀,b_E＝inf(θ(χ))₁＝sup(θ(χ))₁,c_E＝sup(θ(χ))₀，表示模糊随机变量χ的模糊期望值。

②模糊随机水流量被转化成为三角模糊数，系统采用的带乐观-悲观指数的期望值算子将其转化成为确定性值；

首先，定义一个乐观-悲观指数(0≤ε≤1)用于量化决策者的决策偏好，然后使用期望值算子计算该三角模糊数的期望值，并将其作为水流量的确定性变量，则有：

$E^{\mathrm{\ϵ}}\[E^I\[{\widetilde{\stackrel{\‾}{Q}}}_u\]\]=\frac{1-\mathrm{\ϵ}}{2}{a}_E+\frac{1}{2}{b}_E+\frac{\mathrm{\ϵ}}{2}{c}_E;$

其中，表示水流量的确定性变量。

2.河流同化能力模拟。

河流同化能力表示一定的水域、水量条件下，基于维持水域一定的含氧水平，该水域对有机污染物的自净量。其与河段的流量、流速、水温、允许的最低溶解氧含量及污染物本身的特性有关。确定水体对污染物的同化能力是确定地区有机污染物排放标准以及制度合理的污染负荷分配方案的基础。本发明使用Streeter-Phelps模型来模拟河流对有机污染物的同化过程。其主要原理为：当河流受纳有机物后，沿水流方向产生的输移有机物远大于扩散稀释量，当河水流量与污水流量稳定，河水温度不变时，有机物生化降解的好氧量与该时期河水中存在的有机物量成正比。因此，有两个微分耦合模型被定义来描述河流中溶解氧(DO)和生化需氧量(BOD)浓度变化规律。

(1)生化需氧量(BOD)转换机制：生化需氧量表示水中有机物等需氧污染物质含量的一个综合指标。说明水中有机物由于微生物的生化作用进行氧化分解，使之无机化或气体化时所消耗水中溶解氧的总数量。其值越高说明水中有机污染物质越多，污染也就越严重。根据Streeter-Phelps模型的结论，河流中BOD浓度的随时间变化量与其自身浓度成正比，即其中C_b表示BOD浓度，k_d表示脱氧系数。因此，在一个河段t，用k_d(t)表示脱氧系数，d(t)表示该河段长度，v(t)表示该河段平均水流流速，C_b(t)表示流入该河段水流的BOD浓度，基于以上微分模型的解析解，则流出该河段水流的BOD浓度C_b(t+1)为具体如图2所示。定义该浓度等价于：

C_b(t+1)＝α(t)C_b(t)。

(2)溶解氧(DO)转换机制：水中溶解氧的多少是衡量水体自净能力的另一个指标。根据Streeter-Phelps模型的结论，河流中的溶解氧从空气中补充，其补充浓度的随时间变化量与其氧亏(即水体中饱和溶解氧和现存溶解氧的差)浓度成正比，即其中C_b表示BOD浓度，D_o表示氧亏浓度，k_d，k_a表示分别脱氧和充氧系数。与BOD情况相似，在一个河段t，用k_a(t)表示充氧系数，C_b(t)，C_o(t)分别表示流入该河段水流的BOD浓度，C_do表示该河段水流的饱和溶解氧浓度，则流出该河段水流的DO浓度C_o(t+1)为：

$C_o(t+1)=e^{-k_a\left(t\right)\frac{d\left(t\right)}{v\left(t\right)}}{C}_o\left(t\right)-\frac{k_d\left(t\right)}{k_a\left(t\right)-k_d\left(t\right)}\[e^{-k_d\left(t\right)\frac{d\left(t\right)}{v\left(t\right)}}-e^{-k_a\left(t\right)\frac{d\left(t\right)}{v\left(t\right)}}\]C_b\left(t\right)-C_{do}\[1-e^{-k_a\left(t\right)\frac{d\left(t\right)}{v\left(t\right)}}\];$

具体如图2所示。

定义和δ(t)＝C_do[1-β(t)]，该浓度等价于：

C_o(t+1)＝β(t)C_o(t)-γ(t)C_b(t)-δ(t)。

其中，β(t)表示河段t的溶氧浓度的转换系数，γ(t)、δ(t)分别表示影响河段t溶氧浓度的生化需氧量浓度、饱溶氧浓度的转换系数。

3.下层模型构建。

对于下层模型，从河流上游到下游，每个排污控制单元为了达到最小的污染成本，根据地方环保部门制定的污染税收标准以及上游的水质情况依次进行独立的污水BOD削减决策。

(1)河流系统状态方程组构建。

河流系统中的水以及水中的各类化学物质沿着河道顺流而下，在流经过各个的排污控制单元过程中，控制单元取水、排污等活动会一定程度的影响河流状态(水量、各类污染物浓度等)的改变，并且该状态会随水的流动转移到下一个控制单元，如图3所示。因此，在污染负荷分配问题中，需关注河流状态在各个排污之间的转移情况，即流量状态转移方程、BOD浓度状态转移方程、DO浓度状态转移方程被构建。

①流量状态转移方程。基于质量守恒定律以及流量连续性条件，用Q_r(t-1)，Q_r(t)分别表示在监测点t-1和t的流量，用Q_w(t)，Q_d(t)分别表示排污控制单元t的取水量和排污量，则可构建流量状态转移方程如下：

Q_r(t)＝Q_r(t-1)-Q_w(t)+Q_d(t)。

②BOD浓度状态转移方程。基于质量守恒定律、Streeter-Phelps水质模拟模型以及流量状态转移方程，用C_b(t-1)，C_b(t)表示监测点t-1和t的BOD浓度，表示排污控制单元t未处理的污水BOD浓度，以及η(t)为排污控制单元t进行削减污水BOD浓度的比例，其为下层排污者t的决策变量，则可构建BOD浓度状态转移方程如下：

$C_b\left(t\right)Q_r\left(t\right)=\mathrm{\α}(t-1)C_b(t-1)\[Q_r(t-1)-Q_w\left(t\right)\]+C_b^i\left(t\right)\[1-\mathrm{\η}\left(t\right)\]Q_d\left(t\right).$

定义变量则该状态转移方程等价于：

$C_b\left(t\right)=\[1-r\left(t\right)\]\mathrm{\α}(t-1)C_b(t-1)+r\left(t\right)C_b^i\left(t\right)\[1-\mathrm{\η}\left(t\right)\].$

③DO浓度状态转移方程。相似于BOD情形，基于质量守恒定律、Streeter-Phelps水质模拟模型以及流量状态转移方程，用C_o(t-1)，C_o(t)表示监测点t-1和t的DO浓度，表示排污控制单元t未处理的污水DO浓度，则可构建DO浓度状态转移方程如下：

$C_o\left(t\right)Q_r\left(t\right)=\[\mathrm{\β}(t-1)C_o(t-1)-\mathrm{\γ}(t-1)C_b(t-1)-\mathrm{\δ}(t-1)\]\[Q_r(t-1)-Q_w\left(t\right)\]+C_o^i\left(t\right)Q_d\left(t\right)a.$

同样的定义变量则该状态转移方程等价于：

$C_o\left(t\right)=\frac{Q_r(t-1)-Q_w\left(t\right)}{Q_r\left(t\right)}\[\mathrm{\β}(t-1)C_o(t-1)-\mathrm{\γ}(t-1)C_b(t-1)-\mathrm{\δ}(t-1)\]+\frac{Q_d\left(t\right)}{Q_r\left(t\right)}{C}_o^i\left(t\right).$

另外，考虑河流中饱和溶解氧浓度C_do的限制，以上DO浓度状态转移方程可改写为：

$C_o\left(t\right)=min\{\[1-r\left(t\right)\]\[\mathrm{\β}(t-1)C_o(t-1)-\mathrm{\γ}(t-1)C_b(t-1)-\mathrm{\δ}(t-1)\]+r\left(t\right)C_o^i\left(t\right),C_{do}\}.$

(2)初始模型河流状态条件构建。

对于位于整个河流系统中上游第一个排污控制单元，用分别表示上游水流BOD和DO浓度，则该控制单元的模型考虑初始河流状态条件如下：

$\begin{array}{ccc}{Q}_r\left(0\right)=E^{\mathrm{\ϵ}}\[E^I\[{\widetilde{\stackrel{\‾}{Q}}}_u\]\],& C_b\left(0\right)=C_b^u,& C_o\left(0\right)=C_o^u\end{array}.$

(3)下层模型约束条件构建。对于下层模型的决策者-排污控制单元，在进行相应污水BOD浓度削减的决策过程中，需满足一下两个约束条件：允许的BOD和DO浓度标准以及BOD削减能力约束。

①允许的BOD和DO浓度标准。排污控制单元t排放污水到河流中，不能使得在相应检测点河流的BOD浓度(C_b(t))超过环保部门所允许的最大BOD浓度标准(B_u(t))，即：

0≤C_b(t)≤B_u(t)。

同样，不能使相应检测点河流的DO浓度(C_o(t))低于环保部门所允许的最小DO浓度标准(D_l(t))，即：

D_l(t)≤C_o(t)。

②BOD削减能力约束。每个排污控制单元具有污水BOD浓度削减能力的限制，通常采用二级污水处理工艺，如图4所示。对于排污控制单元t，第一级处理后的污水BOD浓度削减率为第二级处理后的污水BOD浓度削减率为在实际的污水处理运行中，环保部门要求所以排污者的污水必须进行第一级处理。因此，排污控制单元t进行污水BOD浓度削减决策，其BOD削减率(η(t))不能够超过同时不能低于即：

${\mathrm{\η}}_t^{\min }\≤\mathrm{\η}\left(t\right)\≤{\mathrm{\η}}_t^{\max }.$

(4)下层模型目标函数构建。下层决策者的目标为最小的单位污染成本目标函数，污染成本由污水处理成本和污染税费组成。

①污水处理成本。基于BOD削减和污水处理成本的经验曲线，如图4所示，本发明对于排污控制单元的污水处理费用采用幂函数形式，用k_c1(t)和k_c2(t)表示排污控制单元t的成本参数1和2，k_ce(t)和k_cb(t)表示排污控制单元t的污水和BOD削减参数，则污水处理成本C(t)可以表示为以下方程：

$C\left(t\right)=k_{c1}\left(t\right)Q_d{\left(t\right)}^{k_{ce}\left(t\right)}+k_{c2}\left(t\right)Q_d{\left(t\right)}^{k_{ce}\left(t\right)}\mathrm{\η}{\left(t\right)}^{k_{cb}\left(t\right)}.$

其次，参数k_c1(t)，k_c2(t)，k_ce(t)和k_cb(t)数值以污水处理构筑物建设和运营所得数据，采用最小二乘法进行曲线拟合所得。

②污染税费。在每个监测点，地方环保部门制定相应的污水税收标准(Al_b(t)，上层模型的决策变量)和税率(T_b)。如果检测点t所对应的河流BOD浓度(C_b(t))低于该税收标准(C_b(t)≤Al_b(t))，则排污控制单元t不需要征收污染税(P(t))；如果检测点t所对应的河流BOD浓度(C_b(t))超过该税收标准(C_b(t)>Al_b(t))，则排污者t需征收污染税P(t)＝Q_r(t)T_b[C_b(t)-Al_b(t)]。因此，排放者t的污染税费模型定义为：

$P\left(t\right)=Q_r\left(t\right)T_b\left\{\frac{|C_b\left(t\right)-{\mathrm{Al}}_b\left(t\right)|+\[C_b\left(t\right)-{\mathrm{Al}}_b\left(t\right)\]}{2}\right\}.$

综上所述，排污控制单元t的目标函数，即最小化的单位污染成本如下：

$\underset{\mathrm{\η}\left(t\right)}{min}{F}_t=\frac{C\left(t\right)}{Q_d\left(t\right)}+\frac{P\left(t\right)}{Q_d\left(t\right)}.$

4.上层模型(地方环保部门)构建。

对于上层模型，地方环保部门为实现公平的流域污染负荷分配，于河流系统中每个监测点制度相应的污染税收标准。

(1)上层模型约束条件构建。对于上层模型的决策者-地方环保部门，在进行相应污染税收标准的制定过程中，需满足允许的BOD浓度标准约束条件。

环保部门所制定的污染税收标准(Al_b(t))不能超过允许的最大BOD浓度标准(B_u(t))，即：

0≤Al_b(t)≤B_u(t)。

(2)上层模型目标函数构建。上层决策者需保证公平的污染负荷分配，因此，本发明采用各排污控制单元的单位污染成本的方差来量化水污染负荷分配的公平性。用表示各排污控制单元平均的单位污染成本，则其最小的方差定义为：

$\underset{{\mathrm{Al}}_b\left(t\right)}{min}dF=\underset{t=1}{\overset{T}{\Σ}}{(F_t-aF)}^2.$

5.全局模型

集成上下层目标函数和上下层约束条件，我们可以得到以下基于污染税的流域污染负荷分配问题的合作型二层多从属者规划全局模型。

$\underset{{\mathrm{Al}}_b\left(t\right)}{min}dF={\mathrm{\Σ}}_{t=1}^T{(F_t-aF)}^2$

如图5所示，表示了针对基于污染税的流域污染负荷分配问题的合作型二层多从属者决策优化建模技术流程图。

针对上述合作型二层多从属者决策模型，本发明提出了混合嵌套基于多智能体的动态极值算法的改进粒子群算法进行求解。由于二层规划问题具有较为特殊的决策规则，即上层决策者只是通过自己的决策去指导下层决策者，并不直接干涉下层的决策；而下层决策者只需要把上层的决策作为参数，他可以在自己的可能范围内自由决策。这种决策机制使得上层决策者在选择策略以优化自己的目标达成时，必须考虑到下层决策者可能采取的策略对自己的不利影响。其求解问题已经为证明为NP-难题，即使是最简单的情况。然而，上述模型中存在非线性的目标函数和约束条件，增加模型计算的难度；另一方面，该模型相比于传统的二层模型具有更为特殊的动态决策结构，即下层模型需要通过河流系统状态转移方程进行依次计算，显然，使得该模型的求解问题变的更加困难。因此，本发明针对基于污染税的流域污染负荷分配二层决策模型，将传统精确算法与现代启发式智能算法相结合，设计了一种混合嵌套粒子群算法，具体包括：①针对上述下层模型特殊的多元动态决策结构，设计一种基于多智能体的动态极值算法为精确快速的获得基于特定上层决策变量的下层模型的最优解；②为计算上层模型，将该算法嵌套至粒子群算法的各个粒子适应值计算的过程中；③为提高整体算法的全局搜索能力，采用Boltzmann选择算子代替传统全局最优粒子的更新机制。下面介绍其主要设计思路。

1.粒子结构设计。为简化描述，以决策排污者为T为例介绍算法的粒子结构设计。在混合嵌套粒子群优化算法中，每个粒子的位置编码为T维的向量，以表示解空间中的一点。每个粒子位子代表上层模型一个决策变量(Al_b(t),t＝1,2,...,T)，如下

$P_1\left(\mathrm{\τ}\right)=\[p_1^1\left(\mathrm{\τ}\right),p_1^2\left(\mathrm{\τ}\right),\mathrm{...},p_1^T\left(\mathrm{\τ}\right)\]\↔\[{\mathrm{Al}}_b\left(1\right),{\mathrm{Al}}_b\left(2\right),\mathrm{...},{\mathrm{Al}}_b\left(T\right)\]$

其中，τ表示迭代次数指标，τ＝1,2,...,I，I表示最大迭代次数，l表示粒子指标，l＝1,2,...,L，L表示种群指标。图8表示了一个包括T个排污控制单元的河流系统的粒子结构表达。

2.粒子性能的混合嵌套计算。基于二层规划特殊的决策规则，在使用粒子群算法进行上层模型的计算的过程中，对各代粒子性能的计算时，需要获得基于特定上层决策变量的下层模型的最优解。因此，针对基于污染税的流域污染负荷分配问题的合作型二层多从属者决策模型中下层模型的多元动态特征，本发明为求得其最优反应函数，设计了一种基于多智能体的动态极值算法。

(1)基于多智能体的动态极值算法。为求解在特定上层决策变量下的上述模型的下层模型，需求解T个一元的非线性规划模型。因此，针对单一的非线性规划模型，本发明提出了一种两阶段极值算法，即第一阶段通过极值算法求解该问题极值，第二阶段通过约束处理方法从第一次获得极值中找到该问题的最优解。在解决单一下层模型求解问题后，本发明基于下层各个模型间的动态传递关系，构建了一种多智能体动态框架，以获得下层T个模型的全部最优解。

①极值算法。以下层第t个单一模型为例，介绍极值算法求解下层单一模型的极值。令下层模型目标函数中的C_b(t)＝Al_b(t)，我们能得到该目标函数的断点为 ${\mathrm{\η}}_p\left(t\right)=1-\frac{{\mathrm{Al}}_b\left(t\right)-\[1-r\left(t\right)\]\mathrm{\α}(t-1)C_b(t-1)}{r\left(t\right)C_b^i\left(t\right)}.$ 因此，该模型的极值存现两种情况：

情况I当C_b(t)≤Al_b(t)时，单一下层模型的目标函数如下：

$F_t^1=k_{c2}\left(t\right)Q_d{\left(t\right)}^{\[k_{ce}\left(t\right)-1\]}\mathrm{\η}{\left(t\right)}^{k_{cb}\left(t\right)}+k_{c1}\left(t\right)Q_d{\left(t\right)}^{\[k_{ce}\left(t\right)-1\]}$

令 $a_t=k_{c2}\left(t\right)Q_d{\left(t\right)}^{\[k_{ce}\left(t\right)-1\]},$ p_t＝k_cb(t)和 $c_t=k_{c1}\left(t\right)Q_d{\left(t\right)}^{\[k_{ce}\left(t\right)-1\]},$ 该目标函数等价于

$F_t^1=a_t\mathrm{\η}{\left(t\right)}^{p_t}+c_t$

因此，在该种情况下，目标函数的极值点为η_e1(t)＝η_p(t)。

情况Ⅱ当C_b(t)>Al_b(t)时，单一下层模型的目标函数如下：

$F_t^2=k_{c2}\left(t\right)Q_d{\left(t\right)}^{\[k_{ce}\left(t\right)-1\]}\mathrm{\η}{\left(t\right)}^{k_{cb}\left(t\right)}-T_b{C}_b^i\left(t\right)\mathrm{\η}\left(t\right)+k_{cl}\left(t\right)Q_d{\left(t\right)}^{\[k_{ce}\left(t\right)-1\]}+T_b\frac{\[1-r\left(t\right)\]\mathrm{\α}(t-1)C_b(t-1)+r\left(t\right)C_b^i\left(t\right)-{\mathrm{Al}}_b\left(t\right)}{r\left(t\right)}.$

令该目标函数等价于：

$F_t^2=a_t\mathrm{\η}{\left(t\right)}^{p_t}-b_t\mathrm{\η}\left(t\right)+c_t+b_t{\mathrm{\η}}_p\left(t\right);$

因此，在该种情况下，该目标函数的极值点为考虑极值点的范围，对该种情况下目标函数的极值点做以下修正：当η_e2(t)≥η_p(t)时，下层模型目标函数的极值点为η_e1(t)＝η_p(t)；当0<η_e2(t)<η_p(t)，下层模型目标函数的极值点为该极值算法的具体步骤如表1所示：

表1：极值算法的步骤：

程序1:求解下层单一模型目标函数极值点
	输入:相关参数；
输出:下层模型目标函数极值点η*(t)；
	步骤1:Letη*(t)＝0；
步骤2:计算模型目标函数断点ηp(t)；
	如果ηp(t)≤0,则到步骤4；

否则到步骤3；
	步骤3:计算情况Ⅱ条件下的断点ηe2(t)；
如果ηe2(t)≥ηp(t),则η*(t)＝ηp(t)；
	否则η*(t)＝ηe2(t)；
步骤4:返回η*(t).

②约束处理方法。在计算下层单一模型目标函数的极值后，需对下层目标进行约束处理，通过对下层模型相应的约束条件的等价转换，允许的BOD浓度标准约束变为：

$\mathrm{\&eta;}\left(t\right)\&GreaterEqual;1-\frac{B_u\left(t\right)-\&lsqb;1-r\left(t\right)\&rsqb;\mathrm{\&alpha;}(t-1)C_b(t-1)}{r\left(t\right)C_b^i\left(t\right)};$

令 ${\mathrm{\&eta;}}_t^{BOD}=1-\frac{B_u\left(t\right)-\&lsqb;1-r\left(t\right)\&rsqb;\mathrm{\&alpha;}(t-1)C_b(t-1)}{r\left(t\right)C_b^i\left(t\right)},$ 该约束可以等价于：

$\mathrm{\&eta;}\left(t\right)\&GreaterEqual;{\mathrm{\&eta;}}_t^{BOD}.$

允许的BOD浓度标准约束变为：

$C_b\left(t\right)\&le;\frac{\mathrm{\&beta;}\left(t\right)C_o\left(t\right)+\mathrm{\&delta;}\left(t\right)-\frac{D_1(t+1)-r(t+1)C_b^i(t+1)}{1-r\left(t\right)}}{r\left(t\right)}.$

令g_t等于上述不等式的右边，该不等式二次转换为：

$\mathrm{\&eta;}\left(t\right)\&GreaterEqual;1-\frac{g_t-\&lsqb;1-r\left(t\right)\&rsqb;\mathrm{\&alpha;}(t-1)C_b(t-1)}{r\left(t\right)C_b^i\left(t\right)}.$

令 ${\mathrm{\&eta;}}_t^{DO}=1-\frac{g_t-\&lsqb;1-r\left(t\right)\&rsqb;\mathrm{\&alpha;}(t-1)C_b(t-1)}{r\left(t\right)C_b^i\left(t\right)},$ 该约束可以等价于：

$\mathrm{\&eta;}\left(t\right)\&GreaterEqual;{\mathrm{\&eta;}}_t^{DO}.$

因此，下层单一模型约束处理方法的具体步骤如表2所示。

表2，下层单一模型约束处理方法的具体步骤：

③多智能体动态框架。基于下层模型之间河流系统动态转移方程，本发明设计了一种基于多智能体的动态框架，结合上述提出的极值算法，以求得在特定上层决策下的全部下层模型的最优反应集。其具体步骤如表3所示。

表3 基于多智能体的动态极值算法的求解过程；

(2)适应值函数。通过上述所提出的算法能够计算获得在特定上层决策变量下的下层模型的最优解，根据双层模型的决策规则，可以得到量化例子性能的适应值函数，即为：

$Fitness\left(P_1\right(\mathrm{\&tau;}\left)\right)=\underset{t=1}{\overset{T}{\&Sigma;}}{(F_t^{{\mathrm{Al}}_b}-{\mathrm{aF}}^{{\mathrm{Al}}_b})}^2;$

其中，因此，整个粒子性能的混合嵌套计算过程可以描述为图7所示的计算方法。

3.粒子更新机制。基于粒子群算法的运行程序，在量化粒子性能后需进行粒子更新，为获得更好性能的粒子。粒子更新方程如下：

P_l(τ+1)＝V_l(τ+1)+P_l(τ)；

其中，表示当前个体最优粒子，即对o＝1,2,...,τ；表示当前粒子群中最优粒子，即对l＝1,2,...,L；θ(τ)表示惯性权重，计算过程为其中θ(1)和θ(I)表示第1和第I代的惯性权重；c_ρ和c_l分别表示个体和全局最优学习因子；r_ρ和r_l表示0到1之间的随机数。为提供粒子的全局搜索能力，本发明采用Boltzmann选择算子代替传统全局最优粒子的更新机制。

(1)Boltzmann选择算子。采用Boltzmann选择算子更新全局最优粒子，其基本思路：每代全局最优粒子通过运用Boltzmann选择算子在各粒子的个体最优粒子中进行选择。为计算各个粒子的个体最优粒子的选择概率，其评价函数定义如下：

${\mathrm{Fit}}_1^B=\frac{1}{1+Fitness\left(P_1^B\right)}.$

其各粒子的个体最优粒子的选择概率为

$\left\{\begin{array}{c}{p}_1\left(\mathrm{\&tau;}\right)=\frac{e^{\left(\frac{{\mathrm{Fit}}_1^B}{{\mathrm{Te}}_{\mathrm{\&tau;}}}\right)}}{\underset{l=1}{\overset{L}{\&Sigma;}}{e}^{\left(\frac{{\mathrm{Fit}}_1^B}{{\mathrm{Te}}_{\mathrm{\&tau;}}}\right)}}\\ T{e}_{\mathrm{\&tau;}}={0.96}^{(\mathrm{\&tau;}-1)}T{e}_0\end{array}\right.$

其中Te₀和Te_τ表示初代和第τ代温度。因此，运用Boltzmann选择算子进行各粒子的全局最优粒子的具体步骤如表4及图8所示。

表4 运用Boltzmann选择算子进行各粒子的全局最优粒子的具体步骤；

4.混合嵌套粒子群算法。该算法的整体步骤如表5所示：

表5 混合嵌套粒子群算法步骤；

为了能进一步了解本发明的发明内容、特点及用法，兹举以如下实例详细说明。

实例：四川省沱江流域某河段为例，其具有5个排污控制单元(i.e.,t＝0,1,2,3,4,5)，整个区域的考虑基于污染税政策下的污染负荷分配决策，该河流系统见图9。通过对该流域水位数据访问调研，该流域饱和溶解氧浓度为10.08mg/L，其他的水位数据见表6。对于各排污控制单元，其污水处理成本参数(k_c1(t)，k_c2(t)，k_ce(t)和k_cb(t))数值以污水处理构筑物建设和运营所得数据，采用最小二乘法进行曲线拟合所得，其他相关数据见表7。对于河流系统的上游情况，其BOD和DO浓度为3.2mg/L和6.2mg/L，其水流量为模糊随机变量，其分布数据见8。假设上层决策的偏好乐观-悲观指数设置为ε＝0.5，通过计算上游水流量的期望值为12.10m³/s。该河流系统各监测点的环境水质标准见表9，另外基本污染税率为0.005yuan·L/mg。

表6 河段水位数据；

表7 排污控单元数据；

表8 河流水流量分别数据；

表9 河流系统各监测点的环境水质标准

计算结果。应用本发明提出的基于污染税的问题的多阶段建模技术对该实例问题进行建模，并采用本发明提出的两步转化方法对模型中的三角模糊随机数进行处理，最后应用本发明提出的混合嵌套粒子群算法对模型进行求解计算，结果如表10所示。通过将算法优化的结果反映，各排污控制单元的污染成本都为1.72yuan·s/m³。因此，地方环保部门通过污染税收的制定，实现了公平的污染负荷分配。这种分配结果可以能够实现流域健康可持续的发展。尽管算法优化的结果是基于模型假设的条件得出的，但在一定程度上仍能反映出实际决策方案尚有提升的空间。在实际决策过程中，各决策者仅仅根据经验判断决定决策方案具有一定风险。运用模型和算法计算的结果可对决策者进行方案选择提供参考，以更好地指导水质管理实践。

表10 模型求解结果

Claims (9)

1.基于污染税的流域污染负荷分配二层决策优化方法，其特征在于，包括以下步骤：

步骤1、系统获取水流量，并将其转换为水流量的确定性变量；

步骤2、系统获取河流中溶解氧浓度关系及生化需氧量浓度关系；

步骤3、系统根据水流量的确定性变量、溶解氧浓度关系及生化需氧量浓度关系构建下层模型；

步骤4、系统根据污染税收标准制定上层模型；

步骤5、系统根据下层模型及上层模型构建基于污染税的流域污染负荷分配二层决策模型；

步骤6、系统对基于污染税的流域污染负荷分配二层决策模型进行求解，得到最优解。

2.根据权利要求1所述的基于污染税的流域污染负荷分配二层决策优化方法，其特征在于，所述步骤1中，首先，系统采用模糊随机变量对未来时间水流量进行量化，系统将水流量分为三个不同的级别:低水位(L)，中水位(M)和高水位(H)，服从概率分布(a,b,c)，三种水位的概率分别可以表示为p_L，p_M和p_H，综上，水流量为:

其次，将水流量转换为模糊随机变量，并计算其模糊期望值，具体包括：

①假设水流量可以表述成为任意一个模糊随机变量χ，含有以下数值：

${\tilde{Q}}_L=(a_L,b_L,c_L),{\tilde{Q}}_M=(a_M,b_M,c_M),{\tilde{Q}}_H=(a_H,b_H,c_H);$

其中，表示三个水流量水平，因此其概率分布如下：

$P\left(\right\{\mathrm{\&omega;}\&Element;\mathrm{\&Omega;}\left|\mathrm{\&chi;}\right(\mathrm{\&omega;})={\tilde{Q}}_L\})=p_L,P\left(\right\{\mathrm{\&omega;}\&Element;\mathrm{\&Omega;}\left|\mathrm{\&chi;}\right(\mathrm{\&omega;})={\tilde{Q}}_M\})=p_M,P\left(\right\{\mathrm{\&omega;}\&Element;\mathrm{\&Omega;}\left|\mathrm{\&chi;}\right(\mathrm{\&omega;})={\tilde{Q}}_H\})=p_H;$

其中，Ω表示三种水流水平及其模糊逻辑；

假设：θ(χ)是模糊随机变量χ的模糊参数，如果θ(χ)＝E^f[X|P],存在，则有：

$inf{\left(\mathrm{\&theta;}\right(\mathrm{\&chi;}\left)\right)}_{\mathrm{\&alpha;}}=E^f\left({\mathrm{inf\&chi;}}_{\mathrm{\&alpha;}}\right|P)=p_1\inf {\left({\tilde{q}}_1\right)}_{\mathrm{\&alpha;}}+p_{2}inf{\left({\tilde{q}}_2\right)}_{\mathrm{\&alpha;}}+p_{3}inf{\left({\tilde{q}}_3\right)}_{\mathrm{\&alpha;}}$

$sup{\left(\mathrm{\&theta;}\right(\mathrm{\&chi;}\left)\right)}_{\mathrm{\&alpha;}}=E^f\left({\mathrm{sup\&chi;}}_{\mathrm{\&alpha;}}\right|P)=p_1\sup {\left({\tilde{q}}_1\right)}_{\mathrm{\&alpha;}}+p_{2}sup{\left({\tilde{q}}_2\right)}_{\mathrm{\&alpha;}}+p_{3}sup{\left({\tilde{q}}_3\right)}_{\mathrm{\&alpha;}};$

那么，模糊随机变量χ的模糊期望值可以表示为则有:

$E^I\&lsqb;\widetilde{\stackrel{\&OverBar;}{Q}}\&rsqb;=(a_E,b_E,c_E);$

其中，a_E＝inf(θ(χ))₀,b_E＝inf(θ(χ))₁＝sup(θ(χ))₁,c_E＝sup(θ(χ))₀，表示模糊随机变量χ的模糊期望值；

②模糊随机水流量被转化成为三角模糊数，系统采用的带乐观-悲观指数的期望值算子将其转化成为确定性值；

首先，定义一个乐观-悲观指数(0≤ε≤1)用于量化决策者的决策偏好，然后使用期望值算子计算该三角模糊数的期望值，并将其作为水流量的确定性变量，则有：

$E^{\mathrm{\&epsiv;}}\&lsqb;E^I\&lsqb;{\widetilde{\stackrel{\&OverBar;}{Q}}}_u\&rsqb;\&rsqb;=\frac{1-\mathrm{\&epsiv;}}{2}{a}_E+\frac{1}{2}{b}_E+\frac{\mathrm{\&epsiv;}}{2}{c}_E;$

其中，表示水流量的确定性变量。

3.根据权利要求1所述的基于污染税的流域污染负荷分配二层决策优化方法，其特征在于，所述步骤2中，系统根据Streeter-Phelps模型的结论，河流中生化需氧量浓度的随时间变化量与其自身浓度成正比，即：

$\frac{{\mathrm{dC}}_b}{dt}=-k_d{C}_b;$

其中，C_b表示生化需氧量浓度，k_d表示脱氧系数；

因此，在一个河段t，用k_d(t)表示脱氧系数，d(t)表示该河段长度，v(t)表示该河段评价水流流速，C_b(t)表示流入该河段水流的生化需氧量浓度，则流出该河段水流的生化需氧量浓度C_b(t+1)为定义该生化需氧量浓度等价于：

C_b(t+1)＝α(t)C_b(t)；

系统根据Streeter-Phelps模型的结论，河流中的溶解氧从空气中补充，其补充浓度的随时间变化量与其氧亏浓度成正比，即：

$\frac{{\mathrm{dD}}_o}{dt}=-k_d{C}_b+k_a{D}_o;$

其中，C_b表示生化需氧量浓度，D_o表示氧亏浓度，k_d表示脱氧系数，k_a表示充氧系数。与生化需氧量情况相似，在一个河段t，用k_a(t)表示充氧系数，C_b(t)，C_o(t)分别表示流入该河段水流的生化需氧量浓度，C_do表示该河段水流的饱和溶解氧浓度，则流出该河段水流的溶解氧浓度C_o(t+1)为：

$C_o(t+1)=e^{-k_a\left(t\right)\frac{d\left(t\right)}{v\left(t\right)}}{C}_o\left(t\right)-\frac{k_d\left(t\right)}{k_a\left(t\right)-k_d\left(t\right)}\&lsqb;e^{-k_d\left(t\right)\frac{d\left(t\right)}{v\left(t\right)}}-e^{-k_a\left(t\right)\frac{d\left(t\right)}{v\left(t\right)}}\&rsqb;C_b\left(t\right)-C_{do}\&lsqb;1-e^{-k_a\left(t\right)\frac{d\left(t\right)}{v\left(t\right)}}\&rsqb;;$

定义： 和δ(t)＝C_do[1-β(t)]，该溶解氧浓度等价于：

C_o(t+1)＝β(t)C_o(t)-γ(t)C_b(t)-δ(t)；

其中，α(t)、β(t)分别可以表示河段t的生化需氧量浓度、溶氧浓度的转换系数，γ(t)、δ(t)分别表示影响河段t溶氧浓度的生化需氧量浓度、饱溶氧浓度转换系数。

4.根据权利要求1所述的基于污染税的流域污染负荷分配二层决策优化方法，其特征在于，所述步骤3中，系统构下层模型包括：

(1)河流系统状态方程组构建，具体包括：

①流量状态转移方程：基于质量守恒定律以及流量连续性条件，用Q_r(t-1)，Q_r(t)分别表示在监测点t-1和t的流量，用Q_w(t)，Q_d(t)分别表示排污控制单元t的取水量和排污量，则可构建流量状态转移方程如下：

Q_r(t)＝Q_r(t-1)-Q_w(t)+Q_d(t)；

②BOD浓度状态转移方程：基于质量守恒定律、Streeter-Phelps水质模拟模型以及流量状态转移方程，用C_b(t-1)，C_b(t)表示监测点t-1和t的BOD浓度，表示排污控制单元t未处理的污水BOD浓度，以及η(t)为排污控制单元t进行削减污水BOD浓度的比例，其为下层排污者t的决策变量，则可构建BOD浓度状态转移方程如下：

$C_b\left(t\right)Q_r\left(t\right)=\mathrm{\&alpha;}(t-1)C_b(t-1)\&lsqb;Q_r(t-1)-Q_w\left(t\right)\&rsqb;+C_b^i\left(t\right)\&lsqb;1-\mathrm{\&eta;}\left(t\right)\&rsqb;Q_d\left(t\right);$

定义变量则该状态转移方程等价于：

$C_b\left(t\right)=\&lsqb;1-r\left(t\right)\&rsqb;\mathrm{\&alpha;}(t-1)C_b(t-1)+r\left(t\right)C_b^i\left(t\right)\&lsqb;1-\mathrm{\&eta;}\left(t\right)\&rsqb;;$

③DO浓度状态转移方程：相似于BOD情形，基于质量守恒定律、Streeter-Phelps水质模拟模型以及流量状态转移方程，用C_o(t-1)，C_o(t)表示监测点t-1和t的DO浓度，表示排污控制单元t未处理的污水DO浓度，则可构建DO浓度状态转移方程如下：

$C_o\left(t\right)Q_r\left(t\right)=\&lsqb;\mathrm{\&beta;}(t-1)C_o(t-1)-\mathrm{\&gamma;}(t-1)C_b(t-1)-\mathrm{\&delta;}(t-1)\&rsqb;\&lsqb;Q_r(t-1)-Q_w\left(t\right)\&rsqb;+C_o^i\left(t\right)Q_d\left(t\right);$

同样的定义变量则该状态转移方程等价于：

$C_o\left(t\right)=\frac{Q_r(t-1)-Q_w\left(t\right)}{Q_r\left(t\right)}\&lsqb;\mathrm{\&beta;}(t-1)C_o(t-1)-\mathrm{\&gamma;}(t-1)C_b(t-1)-\mathrm{\&delta;}(t-1)\&rsqb;+\frac{Q_d\left(t\right)}{Q_r\left(t\right)}{C}_o^i\left(t\right);$

另外，考虑河流中饱和溶解氧浓度C_do的限制，以上DO浓度状态转移方程可改写为：

$C_o\left(t\right)=min\{\&lsqb;1-r\left(t\right)\&rsqb;\&lsqb;\mathrm{\&beta;}(t-1)C_o(t-1)-\mathrm{\&gamma;}(t-1)C_b(t-1)-\mathrm{\&delta;}(t-1)\&rsqb;+r\left(t\right)C_o^i\left(t\right),C_{do}\};$

(2)构建初始模型河流状态条件：

对于位于整个河流系统中上游第一个排污控制单元，用分别表示上游水流BOD和DO浓度，则该控制单元的模型考虑初始河流状态条件如下：

$Q_r\left(0\right)=E^{\mathrm{\&epsiv;}}\&lsqb;E^I\&lsqb;{\widetilde{\stackrel{\&OverBar;}{Q}}}_u\&rsqb;\&rsqb;,C_b\left(0\right)=C_b^u,C_o\left(0\right)=C_o^u;$

(3)构建下层模型约束条件：

①允许的BOD和DO浓度标准：

排污控制单元t排放污水到河流中，不能使得在相应检测点河流的BOD浓度(C_b(t))超过环保部门所允许的最大BOD浓度标准(B_u(t))，即：

0≤C_b(t)≤B_u(t)；

同样，不能使相应检测点河流的DO浓度(C_o(t))低于环保部门所允许的最小DO浓度标准(D_l(t))，即：

D_l(t)≤C_o(t)；

②BOD削减能力约束：

每个排污控制单元具有污水BOD浓度削减能力的限制，通常采用二级污水处理工艺。对于排污控制单元t，第一级处理后的污水BOD浓度削减率为第二级处理后的污水BOD浓度削减率为在实际的污水处理运行中，环保部门要求所以排污者的污水必须进行第一级处理；因此，排污控制单元t进行污水BOD浓度削减决策，其BOD削减率(η(t))不能够超过同时不能低于即：

${\mathrm{\&eta;}}_t^{\min }\&le;\mathrm{\&eta;}\left(t\right)\&le;{\mathrm{\&eta;}}_t^{\max };$

(4)构建下层模型目标函数：

下层决策者的目标为最小的单位污染成本目标函数，污染成本由污水处理成本和污染税费组成；

①污水处理成本：

基于BOD削减和污水处理成本的经验曲线，本发明对于排污控制单元的污水处理费用采用幂函数形式，用k_c1(t)和k_c2(t)表示排污控制单元t的成本参数1和成本参数2，k_ce(t)和k_cb(t)表示排污控制单元t的污水和BOD削减参数，则污水处理成本C(t)可以表示为以下方程：

$C\left(t\right)=k_{c1}\left(t\right)Q_d{\left(t\right)}^{k_{ce}\left(t\right)}+k_{c2}\left(t\right)Q_d{\left(t\right)}^{k_{ce}\left(t\right)}\mathrm{\&eta;}{\left(t\right)}^{k_{cb}\left(t\right)};$

其次，参数k_c1(t)，k_c2(t)，k_ce(t)和k_cb(t)数值以污水处理构筑物建设和运营所得数据，采用最小二乘法进行曲线拟合所得；

②污染税费：

在每个监测点，地方环保部门制定相应的污水税收标准(Al_b(t)，上层模型的决策变量)和税率(T_b)，如果检测点t所对应的河流BOD浓度(C_b(t))低于该税收标准(C_b(t)≤Al_b(t))，则排污控制单元t不需要征收污染税(P(t))；如果检测点t所对应的河流BOD浓度(C_b(t))超过该税收标准(C_b(t)>Al_b(t))，则排污者t需征收污染税P(t)＝Q_r(t)T_b[C_b(t)-Al_b(t)]，因此，排放者t的污染税费模型定义为：

$P\left(t\right)=Q_r\left(t\right)T_b\left\{\frac{|C_b\left(t\right)-A{l}_b\left(t\right)|+\&lsqb;C_b\left(t\right)-A{l}_b\left(t\right)\&rsqb;}{2}\right\};$

综上所述，排污控制单元t的目标函数，即最小化的单位污染成本如下：

$\underset{\mathrm{\&eta;}\left(t\right)}{min}{F}_t=\frac{C\left(t\right)}{Q_d\left(t\right)}+\frac{P\left(t\right)}{Q_d\left(t\right)}.$

5.根据权利要求1所述的基于污染税的流域污染负荷分配二层决策优化方法，其特征在于，所述步骤4中，系统构上层模型包括：

(1)上层模型约束条件构建：

对于上层模型的决策者-地方环保部门，在进行相应污染税收标准的制定过程中，需满足允许的BOD浓度标准约束条件，环保部门所制定的污染税收标准(Al_b(t))不能超过允许的最大BOD浓度标准(B_u(t))，即：

0≤Al_b(t)≤B_u(t)；

(2)上层模型目标函数构建：

上层决策者需保证公平的污染负荷分配，因此，本发明采用各排污控制单元的单位污染成本的方差来量化水污染负荷分配的公平性，用表示各排污控制单元平均的单位污染成本，则其最小的方差定义为：

$\underset{A{l}_b\left(t\right)}{min}dF=\underset{t=1}{\overset{T}{\&Sigma;}}{(F_t-aF)}^2.$

6.根据权利要求1所述的基于污染税的流域污染负荷分配二层决策优化方法，其特征在于，所述步骤5中，系统根据下层模型及上层模型构建的基于污染税的流域污染负荷分配二层决策模型为：

$\underset{A{l}_b\left(t\right)}{min}dF={\mathrm{\&Sigma;}}_{t=1}^T{(F_t-aF)}^2$

7.根据权利要求1所述的基于污染税的流域污染负荷分配二层决策优化方法，其特征在于，所述步骤6中，假设决策排污者为T个，则在特定上层决策变量下地下层模型为T个一元的非线性规划模型；首先，系统利用极值算法计算出单一的非线性规划模型的极值，再通过下层约束条件计算出该单一的非线性规划模型的最优解，进而获得下层T个一元的非线性规划模型的全部最优解集；利用粒子群算法对上述最优解集进行计算，从而获得基于污染税的流域污染负荷分配二层决策模型的最优解。

8.根据权利要求7所述的基于污染税的流域污染负荷分配二层决策优化方法，其特征在于，粒子群算法的粒子更新方程如下：

P_l(τ+1)＝V_l(τ+1)+P_l(τ)；

其中， $P_1^B=\&lsqb;p_1^B\left(1\right),p_1^B\left(2\right),...,p_1^B\left(T\right)\&rsqb;$ 表示当前个体最优粒子，即 $Fitness\left(P_1^B\right)\&le;Fitness\left(P_1\right(o\left)\right),$ 对o＝1,2,...,τ； $P^G=\&lsqb;p_1^G\left(1\right),p_1^G\left(2\right),...,p_1^G\left(T\right)\&rsqb;$ 表示当前粒子群中最优粒子，即 $Fitness\left(P^G\right)\&le;Fitness\left(P_1^B\right),$ 对l＝1,2,...,L；θ(τ)表示惯性权重，计算过程为其中θ(1)和θ(I)表示第1和第I代的惯性权重；c_ρ和分别表示个体和全局最优学习因子；r_ρ和表示0到1之间的随机数。

9.根据权利要求7所述的基于污染税的流域污染负荷分配二层决策优化方法，其特征在于，系统采用Boltzmann选择算子代替传统全局最优粒子的更新机制。