DYNAMIC MULTI-OBJECTIVE PARTICLE SWARM OPTIMIZATION-BASED OPTIMAL CONTROL METHOD FOR WASTEWATER TREATMENT PROCESS

Abstract

A dynamic multi-objective particle swarm optimization based optimal control method is provided to realize the control of dissolved oxygen (S.sub.O) and the nitrate nitrogen (S.sub.NO) in wastewater treatment process. In this method, dynamic multi-objective particle swarm optimization was used to optimize the operation objectives of WWTP, and the optimal solutions of S.sub.O and S.sub.NO can be calculated. Then PID controller was introduced to trace the dynamic optimal solutions of S.sub.O and S.sub.NO. The results demonstrated that the proposed optimal control strategy can address the dynamic optimal control problem, and guarantee the efficient and stable operation. In addition, this proposed optimal control method in this present invention can guarantee the effluent qualities and reduce the energy consumption.

Claims

1. A dynamic multi-objective particle swarm optimization-based optimal control method for a wastewater treatment process (WWTP), comprising: (1) design models of performance indices for the wastewater treatment process: {circle around (1)} analysis dynamic characteristics and operation data of the WWTP, obtain process variables: influent flow rate (Q.sub.in) dissolved oxygen (S.sub.O), nitrate nitrogen (S.sub.NO), ammonia nitrogen (S.sub.NH), suspended solids (SS), which are related to the performance indices including pumping energy (PE), aeration energy (AE) and effluent quality (EQ); {circle around (2)} establish models of the performance indices based on the operation time of S.sub.O and S.sub.NO, the operation time of S.sub.O is thirty minutes, and the operation time of S.sub.NO is two hours, the established models of the performance indices are adjusted per thirty minutes; if an operation time meets the operation time of S.sub.NO only, then the models are designed as: $\begin{array}{cc}\{\begin{array}{c}{f}_1\left(t\right)=\underset{r=1}{\overset{10}{.Math.}}.Math..Math.W_{1.Math..Math.r}\left(t\right)e^{-{.Math.x\left(t\right)-c_{1.Math..Math.r}\left(t\right).Math.}^2/2.Math.{b_{1.Math..Math.r}\left(t\right)}^2}+W_1\left(t\right)\\ f_2\left(t\right)=\underset{r=1}{\overset{10}{.Math.}}.Math..Math.W_{2.Math..Math.r}\left(t\right)e^{-{.Math.x\left(t\right)-c_{2.Math..Math.r}\left(t\right).Math.}^2/2.Math.{b_{2.Math..Math.r}\left(t\right)}^2}+W_2\left(t\right)\end{array}& \left(1\right)\end{array}$ where f.sub.1(t) is PE model at tth time, f.sub.2(t) is EQ model at tth time, e.sup.x(t)c.sup.1r.sup.(t).sup.2.sup./2b.sup.1r.sup.(t).sup.2 and e.sup.1(t)c.sup.2r.sup.(t).sup.2.sup./2b.sup.2r.sup.(t).sup.2 are rth radial basis function of f.sub.1(t) and f.sub.2(t) at the tth time, r=1, 2, . . . , 10, x(t)=[Q.sub.in(t), S.sub.O(t), S.sub.NO(t), S.sub.NH(t), SS(t)] is an input vector of PE model and EQ model, c.sub.1r(t) and c.sub.2r(t) are centers of the rth radial basis function of f.sub.1(t) and f.sub.2(t) at the tth time, and ranges of c.sub.1r(t) and c.sub.2r(t) are [1, 1] respectively; b.sub.1r(t) and b.sub.2r(t) are widths of the rth radial basis function of f.sub.1(t) and f.sub.2(t) at the tth time, and ranges of b.sub.1r(t) and b.sub.2r(t) are [0, 2] respectively, W.sub.1r(t) and W.sub.2r(t) are weights of the rth radial basis function of f.sub.1(t) and f.sub.2(t) at the tth time, and ranges of W.sub.1r(t) and W.sub.2r(t) are [3, 3] respectively; W.sub.1(t) and W.sub.2(t) are output offsets of the rth radial basis function of f.sub.1(t) and f.sub.2(t) at the tth time, and ranges of W.sub.1(t) and W.sub.2(t) are [2, 2] respectively; if the operation time meets the time of S.sub.O, then the models are designed as: $\begin{array}{cc}\{\begin{array}{c}{f}_1\left(t\right)=\underset{r=1}{\overset{10}{.Math.}}.Math..Math.W_{1.Math..Math.r}\left(t\right)e^{-{.Math.x\left(t\right)-c_{1.Math..Math.r}\left(t\right).Math.}^2/2.Math.{b_{1.Math..Math.r}\left(t\right)}^2}+W_1\left(t\right)\\ f_2\left(t\right)=\underset{r=1}{\overset{10}{.Math.}}.Math..Math.W_{2.Math..Math.r}\left(t\right)e^{-{.Math.x\left(t\right)-c_{2.Math..Math.r}\left(t\right).Math.}^2/2.Math.{b_{2.Math..Math.r}\left(t\right)}^2}+W_2\left(t\right)\\ {}_{f_3}.Math.\left(t\right)=\underset{r=1}{\overset{10}{.Math.}}.Math..Math.W_{3.Math..Math.r}\left(t\right)e^{-{.Math.x\left(t\right)-c_{3.Math..Math.r}\left(t\right).Math.}^2/2.Math.{b_{3.Math..Math.r}\left(t\right)}^2}+W_3\left(t\right)\end{array}& \left(2\right)\end{array}$ where f.sub.3(t) is AE model at the tth time, e.sup.s(t)c.sup.3r.sup.(t).sup.2.sup./2b.sup.3r.sup.(t).sup.2 is rth radial basis function of f.sub.3(t) at the tth time, c.sub.3r(t) is center of the rth radial basis function of f.sub.3(t) at the tth time, and a range of c.sub.3r(t) is [1, 1]; b.sub.3r(t) is widths of the rth radial basis function of f.sub.3(t) at the tth time, and a range of b.sub.3r(t) is [0, 2]; W.sub.3r(t) is weights of the rth radial basis function of f.sub.3(t) at the tth time, and a range of W.sub.3r(t) is [3, 3]; W.sub.3(t) is an output offset of the rth radial basis function of f.sub.3(t), and a range of W.sub.3(t) is [2, 2]; (2) dynamic optimization of the control variables for WWTP (2)-1 set maximum iterative numbers of optimization process T.sub.max; (2)-2 take the established models of the performance indices as optimization objectives; (2)-3 regard inputs of the optimization objectives x(t)=[Q.sub.in(t), S.sub.O(t), S.sub.NO(t), S.sub.NH(t), SS(t)] as the position of particles, calculate values of the optimization objectives, update personal optimal position (pBest.sub.k,i(t) and the position and velocity of the particles, the update process is:
x.sub.k,i(t+1)=x.sub.k,i(t)+v.sub.k,i(t+1) (3)
v.sub.k,i(t+1)=(t)v.sub.k,i(t)+c.sub.1.sub.1(pBest.sub.k,i(t)x.sub.k,i(t))+c.sub.2.sub.2(gBest.sub.k(t)(t)) (4) where x.sub.k,i(t+1) is the position of ith particle in kth iteration at t+1th time, v.sub.k,i(t+1) is the velocity of the ith particle in the kth iteration at the t+1th time, is inertia weight, a range of is [0, 1], cl.sub.1 and c.sub.2 are learning parameters, .sub.1 and .sub.2 are uniformly distributed random numbers, pBest.sub.k,i(t) is personal optimal position of the ith particle in the kth iteration at the tth time, and gBest.sub.k(t) is personal optimal position in the kth iteration at the tth time; (2)-4 design diversity index and convergence index based on the Chebyshev distance, the diversity index is designed to measure distribution quality of non-dominated solutions, $\begin{array}{cc}U\left(t\right)=\sqrt{\frac{1}{\mathrm{NS}\left(t\right)-1}.Math.\underset{m=1}{\overset{\mathrm{NS}\left(t\right)}{.Math.}}.Math..Math.{\left(\stackrel{.Math.}{D}\left(t\right)-D_m\left(t\right)\right)}^2}& \left(5\right)\end{array}$ where U(t) is diversity of optimal solutions at the tth time, m=1, 2, . . . , NS(t), NS(t) is the number of non-dominated solutions at time t, (t) is average distance of all the Chebyshev distance D.sub.m(t), D.sub.m(t) is Chebyshev distance of consecutive solutions of the mth solution; and the convergence index is developed to obtain degree of proximity, which is calculated as $\begin{array}{cc}A\left(t\right)=\frac{1}{\mathrm{NS}\left(t\right)}.Math.\sqrt{\underset{l=1}{\overset{\mathrm{NS}\left(t\right)}{.Math.}}.Math..Math.d_l\left(t\right)}& \left(6\right)\end{array}$ where A(t) is the convergence of the optimal solutions at the tth time, d.sub.l(t) is the Chebyshev distance of the lth solution between the kth iteration and the k1th iteration; (2)-5 judge changes of the optimization objectives, if the number of the objectives is changed, return to step (2)-6; otherwise, return to (2)-7; (2)-6 when the number of the objectives is increased, some particles will be changed to enhance diversity performance, the update process of population size is $\begin{array}{cc}{N}_{k+1}\left(t\right)=\{\begin{array}{cc}{N}_k\left(t\right)& {}_k\left(t\right)=0\\ N_k\left(t\right)-\left(N_k\left(t\right)-{\mathrm{NS}}_k\left(t\right)\right).Math.{}_k\left(t\right)& {}_k\left(t\right)<0\\ N_k\left(t\right)+{\mathrm{NS}}_k\left(t\right).Math.{}_k\left(t\right)& {}_k\left(t\right)>0\end{array}& \left(7\right)\end{array}$ where N.sub.k+1(t) and N.sub.k(t) are population size in kth iteration and in k+1th iteration at the tth time respectively, .sub.k(t) is gradient of diversity in the kth iteration at the tth time, which is calculated as $\begin{array}{cc}{}_k\left(t\right)=\frac{U_k\left(t\right)-U_k\left(t-\mathrm{.Math.}\right)}{\mathrm{.Math.}}& \left(8\right)\end{array}$ where is an adjusted frequency of the population size, and a range of is [1, T.sub.max]; if the number of the objectives is decreased, some particles will be changed to improve the convergence performance, the update process of the population size is $\begin{array}{cc}{N}_{k+1}\left(t\right)=\{\begin{array}{cc}{N}_k\left(t\right)& {}_k\left(t\right)=0\\ N_k\left(t\right)+{\mathrm{NS}}_k\left(t\right).Math.{}_k\left(t\right)& {}_k\left(t\right)<0\\ N_k\left(t\right)-\left(N_k\left(t\right)-{\mathrm{NS}}_k\left(t\right)\right).Math.{}_k\left(t\right)& {}_k\left(t\right)>0\end{array}& \left(9\right)\end{array}$ where .sub.k(t) is gradient of convergence in the kth iteration at the tth time, which is calculated as $\begin{array}{cc}{}_k\left(t\right)=\frac{A_k\left(t\right)-A_k\left(t-\mathrm{.Math.}\right)}{\mathrm{.Math.}}& \left(10\right)\end{array}$ (2)-7 compare pBest.sub.k(t) with solutions .sub.k1(t) in the archive, where .sub.k1(t)=[.sub.k1,1(t), .sub.k1,2(t) , . . . , .sub.k1,.Math.(t)], .sub.k1,.Math.(t) is .Math.th optimal solutions in k1th iteration at the tth time of the archive, the archive .sub.k(t) is updated by a dominated relationship, and the calculation process of the dominated relationship is:
.sub.k(t)=.sub.k1(t)p.sub.k1(t), if f.sub.h(a.sub.k.Math.(t))f.sub.h(p.sub.k(t)), h=1,2,3 (11) where is a relationship of combine, if a value of pBest.sub.k1(t) is lower than an objective value of a.sub.k1,.Math.(t), then the pBest.sub.k1(t) will be saved in the archive, otherwise, a.sub.k1,.Math.(t) will be saved, then gBest.sub.k(t) will be selected from the archive according to the density method; (2)-8 if the current iteration is greater than the preset T.sub.max, then return to step (2)-9, otherwise, return to step (2)-3; (2)-9 select a set of global optimal solutions gBest.sub.Tmax(t) from the archive randomly, and gBest.sub.Tmax(t)=[Q.sub.in,Tmax*(t), S.sub.O,Tmax*(t), S.sub.NO,Tmax*(t), S.sub.NH,Tmax*(t), SS.sub.Tmax*(t)], where Q.sub.in,Tmax*(t) is an optimal solution of influent flow rate, S.sub.O,Tmax*(t) is an optimal solution of dissolved oxygen, S.sub.NO,Tmax*(t) is an optimal solution of nitrate nitrogen, S.sub.NH,Tmax*(t) is an optimal solution of ammonia nitrogen, SS.sub.Tmax*(t) is an optimal solution of suspend solid; (3) tracking control of the optimal solutions in WWTP (3)-1 design the multivariable PID controller, the output of PID controller is shown as $\begin{array}{cc}.Math..Math.u\left(t\right)=K_p\left[e\left(t\right)+H_{}.Math.{}_0^t.Math.e\left(t\right).Math.\mathrm{dt}+H_d.Math.\frac{\mathrm{de}\left(t\right)}{\mathrm{dt}}\right]& \left(12\right)\end{array}$ where u(t)=[K.sub.La.sub.5(t), Q.sub.a(t)].sup.T, K.sub.La.sub.5(t) is error of oxygen transfer coefficient in fifth unit at time t, Q.sub.a(t) is error of internal recycle flow rate at time t, K.sub.p is a proportionality coefficient; H.sub. is a integral time constant; H.sub.d is a differential time constant; e(t) is error between a real output and the optimal solution
e(t)=z(t)y(t) (13) where e(t)=[e.sub.1(t), e.sub.2(t)].sup.T, e.sub.1(t) and e.sub.2(t) are errors of S.sub.O and S.sub.NO, respectively; z(t)=[z.sub.1(t), z.sub.2(t)].sup.T, z.sub.1(t) is an optimal set-point concentration of S.sub.O at time t, z.sub.2(t) is an optimal set-point concentration of S.sub.NO at time t, y(t)=[y.sub.1(t), y.sub.2(t)].sup.T, y.sub.1(t) is the concentration of S.sub.O at time t, y.sub.2(t) is the concentration of S.sub.NO at time t; (3)-2 outputs of PID controller are variation of manipulated variables oxygen transfer coefficient (K.sub.La) and internal circulation return flow (Q.sub.a); (4) take K.sub.La and Q.sub.a as an input of a control system of WWTP, and then control S.sub.O and S.sub.NO by the calculated K.sub.La and Q.sub.a, and outputs of the control system in WWTP are real concentrations of S.sub.O and S.sub.NO.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0028] FIG. 1 shows the results of S.sub.O in DMOPSP-based method.

[0029] FIG. 2 shows the results of S.sub.NO in DMOPSP-based method.

DETAILED DESCRIPTION

[0030] A dynamic multi-objective particle swarm optimization-based optimal control method is designed for wastewater treatment process (WWTP), the steps are:

[0031] (1) design the models of performance indices for wastewater treatment process:

[0032] {circle around (1)} analysis the dynamic characteristics and the operation data of WWTP, obtain the related process variables of the key performance indices pumping energy (PE), aeration energy (AE) and effluent quality (EQ): influent flow rate (Q.sub.in), dissolved oxygen (S.sub.O), nitrate nitrogen (S.sub.NO), ammonia nitrogen (S.sub.NH), suspended solids (SS);

[0033] {circle around (2)} establish the models of the performance indices based on the operation time of S.sub.O and S.sub.NO, the operation time of S.sub.O is thirty minutes, and the operation time of S.sub.NO is two hours, the established models of the performance indices are adjusted per thirty minutes; if the operation time meets the operation time of S.sub.NO, then the models are designed as:

[00010]$\begin{array}{cc}\{\begin{array}{c}{f}_1\left(t\right)=\underset{r=1}{\overset{10}{.Math.}}.Math..Math.W_{1.Math..Math.r}\left(t\right)e^{-{.Math.x\left(t\right)-c_{1.Math..Math.r}\left(t\right).Math.}^2/2.Math.{b_{1.Math..Math.r}\left(t\right)}^2}+W_1\left(t\right)\\ f_2\left(t\right)=\underset{r=1}{\overset{10}{.Math.}}.Math..Math.W_{2.Math..Math.r}\left(t\right)e^{-{.Math.x\left(t\right)-c_{2.Math..Math.r}\left(t\right).Math.}^2/2.Math.{b_{2.Math..Math.r}\left(t\right)}^2}+W_2\left(t\right)\end{array}& \left(1\right)\end{array}$

where f.sub.1(t) is PE model at the tth time, f.sub.2(t) is EQ model at the tth time, e.sup.x(t)c.sup.1r.sup.(t).sup.2.sup./2b.sup.1r.sup.(t).sup.2 and e.sup.1(t)c.sup.2r.sup.(t).sup.2.sup./2b.sup.2r.sup.(t).sup.2 are the rth radial basis function of f.sub.1(t) and f.sub.2(t) at the tth time, r=1, 2, . . . , 10, x(t)=[Q.sub.in(t), S.sub.O(t), S.sub.NO(t), S.sub.NH(t), SS(t)] is the input vector of PE model and EQ model, c.sub.1r(t) and c.sub.2r(t) are the centers of the rth radial basis function of f.sub.1(t) and f.sub.2(t) at the tth time, and the ranges of c.sub.1r(t) and c.sub.2r(t) are [1, 1] respectively; b.sub.1r(t) and b.sub.2r(t) are the widths of the rth radial basis function of f.sub.1(t) and f.sub.2(t) at the tth time, and b.sub.1r(t)=1.2, b.sub.2r(t)=1.2; W.sub.1r(t) and W.sub.2r(t) are the weights of the rth radial basis function of f.sub.1(t) and f.sub.2(t) at the tth time, and W.sub.1r(t)=2.4, W.sub.2r(t)=1.8; W.sub.1(t) and W.sub.2(t) are the output offsets of the rth radial basis function of f.sub.1(t) and f.sub.2(t) at the tth time, and W.sub.1(t)=0.8, W.sub.2(t)=0.2; if the operation time meets the time of S.sub.O, then the models are designed as:

[00011]$\begin{array}{cc}\{\begin{array}{c}{f}_1\left(t\right)=\underset{r=1}{\overset{10}{.Math.}}.Math..Math.W_{1.Math..Math.r}\left(t\right)e^{-{.Math.x\left(t\right)-c_{1.Math..Math.r}\left(t\right).Math.}^2/2.Math.{b_{1.Math..Math.r}\left(t\right)}^2}+W_1\left(t\right)\\ f_2\left(t\right)=\underset{r=1}{\overset{10}{.Math.}}.Math..Math.W_{2.Math.r}\left(t\right)e^{-{.Math.x\left(t\right)-c_{2.Math..Math.r}\left(t\right).Math.}^2/2.Math.{b_{2.Math..Math.r}\left(t\right)}^2}+W_2\left(t\right)\\ f_3\left(t\right)=\underset{r=1}{\overset{10}{.Math.}}.Math..Math.W_{3.Math..Math.r}\left(t\right)e^{-{.Math.x\left(t\right)-c_{3.Math..Math.r}\left(t\right).Math.}^2/2.Math.{b_{3.Math.r}\left(t\right)}^2}+W_3\left(t\right)\end{array}& \left(2\right)\end{array}$

where f.sub.3(t) is AE model at the tth time, e.sup.s(t)c.sup.3r.sup.(t).sup.2.sup./2b.sup.3r.sup.(t).sup.2 is the rth radial basis function of f.sub.3(t) at the tth time, c.sub.3r(t) is the center of the rth radial basis function of f.sub.3(t) at the tth time, and the range of c.sub.3r(t) is [1, 1]; b.sub.3r(t) is the widths of the rth radial basis function of f.sub.3(t) at the tth time, and b.sub.3r(t)=0.5, W.sub.3r(t) is the weights of the rth radial basis function of f.sub.3(t) at the tth time, and W.sub.3r=1.4; W.sub.3(t) is the output offset of the rth radial basis function of f.sub.3(t), and W.sub.2(t)=1.2;

[0034] (2) dynamic optimization of the control variables for WWTP

[0035] (2)-1 set the maximum iterative numbers of the optimization process T.sub.max, T.sub.max=200;

[0036] (2)-2 take the established models of performance indices as optimization objectives;

[0037] (2)-3 regard the inputs of the optimization objectives x(t)=[Q.sub.in(t), S.sub.O(t), S.sub.NO(t), S.sub.NH(t), SS(t)] as the position of the particles, calculate values of the optimization objectives, update the personal optimal position (pBest.sub.k,i(t)) and the position and the velocity of the particles, the update process are:

x.sub.k,i(t+1)=x.sub.k,i(t)+v.sub.k,i(t+1) (3)

v.sub.k,i(t+1)=(t)v.sub.k,i(t)+c.sub.1.sub.1(pBest.sub.k,i(t)x.sub.k,i(t))+c.sub.2.sub.2(gBest.sub.k(t)x.sub.k,i(t)) (4)

where x.sub.k,i(t+1) is the position of the ith particle in the kth iteration at the t+1th time, v.sub.k,i(t+1) is the velocity of the ith particle in the kth iteration at the t+1th time, is the inertia weight, =0.8, c.sub.1 and c.sub.2 are the learning parameters, c.sub.1=0.4, c.sub.2=0.4, .sub.1 and .sub.2 are the uniformly distributed random numbers, .sub.1=0.2, .sub.2=0.2, pBest.sub.k,i(t) is the personal optimal position of the ith particle in the kth iteration at the tth time, and gBest.sub.k(t) is the personal optimal position in the kth iteration at the tth time;

[0038] (2)-4 design the diversity index and the convergence index based on the Chebyshev distance, diversity index is designed to measure the distribution quality of the non-dominated solutions,

[00012]$\begin{array}{cc}U\left(t\right)=\sqrt{\frac{1}{\mathrm{NS}\left(t\right)-1}.Math.\underset{m=1}{\overset{\mathrm{NS}\left(t\right)}{.Math.}}.Math..Math.{\left(\stackrel{.Math.}{D}\left(t\right)-D_m\left(t\right)\right)}^2}& \left(5\right)\end{array}$

where U(t) is the diversity of the optimal solutions at the tth time, m=1, 2, . . . , NS(t), NS(t) is the number of non-dominated solutions at time t, (t) is the average distance of all the Chebyshev distance D.sub.m(t), D.sub.m(t) is the Chebyshev distance of consecutive solutions of the mth solution; and convergence index is developed to obtain the degree of proximity, which are calculated as

[00013]$\begin{array}{cc}A\left(t\right)=\frac{1}{\mathrm{NS}\left(t\right)}.Math.\sqrt{\underset{l=1}{\overset{\mathrm{NS}\left(t\right)}{.Math.}}.Math..Math.d_l\left(t\right)}& \left(6\right)\end{array}$

where A(t) is the convergence of the optimal solutions at the tth time, d.sub.l(t) is the

[0039] Chebyshev distance of the lth solution between the kth iteration and the k1th iteration;

[0040] (2)-5 judge the changes of the optimization objectives, if the number of the objectives is changed, return to step (2)-6; otherwise, return to (2)-7;

[0041] (2)-6 when the number of the objectives is increased, some particles will be changed to enhance the diversity performance, the update process of the population size is

[00014]$\begin{array}{cc}{N}_{k+1}\left(t\right)=\{\begin{array}{cc}{N}_k\left(t\right)& {}_k\left(t\right)=0\\ N_k\left(t\right)-\left(N_k\left(t\right)-{\mathrm{NS}}_k\left(t\right)\right).Math.{}_k\left(t\right)& {}_k\left(t\right)<0\\ N_k\left(t\right)+{\mathrm{NS}}_k\left(t\right).Math..Math.{}_k\left(t\right)& {}_k\left(t\right)>0\end{array}& \left(7\right)\end{array}$

where N.sub.k+1(t) and N.sub.k(t) are the population size in the kth iteration and in the k+1th iteration at the tth time respectively, .sub.k(t) is the gradient of diversity in the kth iteration at the tth time, which is calculated as

[00015]$\begin{array}{cc}{}_k\left(t\right)=\frac{U_k\left(t\right)-U_k\left(t-\mathrm{.Math.}\right)}{\mathrm{.Math.}}& \left(8\right)\end{array}$

where is the adjusted frequency of population size, and the range of is [1, 200]; if the number of the objectives is decreased, some particles will be changed to improve the convergence performance, the update process of the population size is

[00016]$\begin{array}{cc}{N}_{k+1}\left(t\right)=\{\begin{array}{cc}{N}_k\left(t\right)& {}_k\left(t\right)=0\\ N_k\left(t\right)+{\mathrm{NS}}_k\left(t\right).Math.{}_k\left(t\right)& {}_k\left(t\right)<0\\ N_k\left(t\right)-\left(N_k\left(t\right)-{\mathrm{NS}}_k\left(t\right)\right).Math.{}_k\left(t\right)& {}_k\left(t\right)>0\end{array}& \left(9\right)\end{array}$

where .sub.k(t) is the gradient of convergence in the kth iteration at the tth time, which is calculated as

[00017]$\begin{array}{cc}{}_k\left(t\right)=\frac{A_k\left(t\right)-A_k\left(t-\mathrm{.Math.}\right)}{\mathrm{.Math.}}& \left(10\right)\end{array}$

[0042] (2)-7 compare pBest.sub.k(t) with the solutions .sub.k1(t) in the archive, where .sub.k1(t)=[.sub.k1,1(t), .sub.k1,2(t), . . . , .sub.k1,.Math.(t)], .sub.k1,.Math.(t) is the .Math.th optimal solutions in the k-1th iteration at the tth time of the archive, the archive .sub.k(t) is updated by the dominated relationship, and the calculation process of the dominated relationship can be shown as

.sub.k(t)=.sub.k1(t)p.sub.k1(t), if f.sub.h(a.sub.ki(t))f.sub.h(p.sub.k(t)), h=1, 2, 3 (11)

where is the relationship of combine, if the value of pBest.sub.k-1(t) is lower than the objective value of a.sub.k1,.Math.(t), then the pBest.sub.k1(t) will be saved in the archive, otherwise, a.sub.k1,i.Math.(t) will be saved, then gBest.sub.k(t) will be selected from the archive according to the density method;

[0043] (2)-8 if the current iteration is greater than the preset T.sub.max, then return to step (2)-9, otherwise, return to step (2)-3;

[0044] (2)-9 select a set of global optimal solutions gBest.sub.Tmax(t) from the archive randomly, and gBest.sub.Tmax(t)=[Q.sub.in,Tmax*(t), S.sub.O,Tmax*(t), S.sub.NO,Tmax*(t), S.sub.NH,Tmax*(t), SS.sub.Tmax*(t)], where Q.sub.in,Tmax* (t) is the optimal solution of influent flow rate, S.sub.O,Tmax*(t) is the optimal solution of dissolved oxygen, S.sub.NO,Tmax*(t) is the optimal solution of nitrate nitrogen, S.sub.NH,Tmax*(t) is the optimal solution of ammonia nitrogen, SS.sub.Tmax*(t) is the optimal solution of suspend solid;

[0045] (3) tracking control of the optimal solutions in WWTP

[0046] (3)-1 design the multivariable PID controller, the output of PID controller is shown as

[00018]$\begin{array}{cc}.Math..Math.u\left(t\right)=K_p\left[e\left(t\right)+H_{}.Math.{}_0^t.Math.e\left(t\right).Math.\mathrm{dt}+H_d.Math.\frac{\mathrm{de}\left(t\right)}{\mathrm{dt}}\right]& \left(12\right)\end{array}$

where u(t)=[K.sub.La.sub.5(t), Q.sub.a(t)].sup.T, K.sub.La.sub.5(t) is the error of the oxygen transfer coefficient in the fifth unit at time t, Q.sub.a(t) is the error of the internal recycle flow rate at time t, K.sub.p is the proportionality coefficient; H.sub. is the integral time constants; H.sub.d is the differential time constants; e(t) is the error between the real output and the optimal solution

e(t)=z(t)y(t) (13)

where e(t)=[e.sub.1(t), e.sub.2(t)].sup.T, e.sub.1(t) and e.sub.2(t) are the errors of S.sub.O and S.sub.NO, respectively; z(t)=[z.sub.1(t), z.sub.2(t)].sup.T, z.sub.1(t) is the optimal set-point concentration of S.sub.O at time t, z.sub.2(t) is the optimal set-point concentration of S.sub.NO at time t, y(t)=[y.sub.1(t), y.sub.2(t)].sup.T, y.sub.1(t) is the concentration of S.sub.O at time t, y.sub.2(t) is the concentration of S.sub.NO at time t;

[0047] (3)-2 the outputs of PID controller are the variation of the manipulated variables oxygen transfer coefficient (K.sub.La) and internal circulation return flow (Q.sub.a);

[0048] (4) take K.sub.La and Q.sub.a as the input of the control system of WWTP, and then control S.sub.O and S.sub.NO by the calculated K.sub.La and Q.sub.a, and the outputs of the control system in WWTP are the real concentration of S.sub.O and S.sub.NO.

[0049] The outputs of the control system based on DMOPSO-based optimal control method is the concentration of S.sub.O and S.sub.NO. FIG. 1(a) gives S.sub.O values, X axis shows the time, and the unit is day, Y axis is control results of S.sub.O, and the unit is mg/L. FIG. 1(b) gives the control errors of S.sub.O, X axis shows the time, and the unit is day, Y axis is control errors of S.sub.O, and the unit is mg/L. FIG. 2(a) gives the S.sub.NO values, X axis shows the time, and the unit is day, Y axis is control results of S.sub.NO, and the unit is mg/L. FIG. 2(b) gives the control errors of S.sub.NO, X axis shows the time, and the unit is day, Y axis is control errors of S.sub.NO, and the unit is mg/L.