COOPERATIVE OPTIMAL CONTROL METHOD AND SYSTEM FOR WASTEWATER TREATMENT PROCESS

Abstract

In a cooperative optimal control system, firstly, two-level models are established to capture the dynamic features of different time-scale performance indices. Secondly, a data-driven assisted model based cooperative optimization algorithm is developed to optimize the two-level models, so that the optimal set-points of dissolved oxygen and nitrate nitrogen can be acquired. Thirdly, a predictive control strategy is designed to trace the obtained optimal set-points of dissolved oxygen and nitrate nitrogen. This proposed cooperative optimal control system can effectively deal with the difficulties of formulating the dynamic features and acquiring the optimal set-points.

Claims

1. A cooperative optimal control method for wastewater treatment process (WWTP), comprising the following steps: (1) select process variables of pumping energy, which include nitrate nitrogen (S.sub.NO), mixed liquor suspended solids (MLSS), and choose process variables of aeration energy and the effluent quality, which include dissolved oxygen (S.sub.O), suspended solids (SS), ammonia nitrogen (S.sub.NH), S.sub.NO; (2) formulate two-level models based on different time scales, wherein an upper-level model is designed for the pumping energy, and lower-level models are designed for the aeration energy and effluent quality:
F.sub.1(t.sub.1)=l.sub.1(x.sub.u(t.sub.1)),(1)
f.sub.1(t.sub.2)=l.sub.2(x.sub.l(t.sub.2),x.sub.u*(t.sub.1))
f.sub.2(t.sub.2)=l.sub.3(x.sub.l(t.sub.2),x.sub.u*(t.sub.1)),(2) where F.sub.1(t.sub.1) is a model of the pumping energy at time t.sub.1, l.sub.1(x.sub.u(t.sub.1)) is a mapping function of the pumping energy model, f.sub.1(t.sub.2) is a model of the aeration energy at time t.sub.2, l.sub.2(x.sub.l(t.sub.2), x*.sub.u(t.sub.1)) is a mapping function of the aeration energy model, x*.sub.u(t.sub.1) is an optimal set-points of nitrate nitrogen S.sub.NO* at time t.sub.1, f.sub.2(t.sub.2) is a model of the effluent quality at time t.sub.2, l.sub.3(x.sub.l(t.sub.2), x*.sub.u(t.sub.1)) is a mapping function of the effluent quality model, x.sub.u(t.sub.1)=[S.sub.NO(t.sub.1), MLSS(t.sub.1)] is input variables vector of the pumping energy at time t.sub.1, S.sub.NO(t.sub.1) is concentration of nitrate nitrogen S.sub.NO at time t.sub.1, MLSS(t.sub.1) is concentration of the mixed liquor suspended solids MLSS at time t.sub.1, x.sub.l(t.sub.2)=[S.sub.O(t.sub.2), SS(t.sub.2), S.sub.NH(t.sub.2)], S.sub.O(t.sub.2) is concentration of S.sub.O at time t.sub.2, SS(t.sub.2) is concentration of SS at time t.sub.2, S.sub.NH(t.sub.2) is concentration of S.sub.NH at time t.sub.2, and [S.sub.O(t.sub.2), SS(t.sub.2), S.sub.NH(t.sub.2), S.sub.NO*(t.sub.1)] are input variables vector of the aeration energy and effluent quality at time t.sub.2; (3) design a cooperative optimization algorithm to optimize upper-level and lower-level optimization problems to obtain optimal set-points of control variables, an optimization period of upper level is two hours, and an optimization period of lower level is thirty minutes, which includes: {circle around (1)} formulate the upper-level and lower-level problems:
Min F.sub.1(S.sub.NO(t.sub.1),MLSS(t.sub.1)),(3)
Min[f.sub.1(S.sub.O(t.sub.2),S.sub.NH(t.sub.2),SS(t.sub.2),S.sub.NO*(t.sub.1)),
f.sub.2(S.sub.O(t.sub.2),S.sub.NH(t.sub.2),SS(t.sub.2),S.sub.NO*(t.sub.1))],(4) where Min F.sub.1(S.sub.NO(t.sub.1), MLSS(t.sub.1)) is the upper-level optimization problem, Min [f.sub.1(S.sub.O(t.sub.2), S.sub.NH(t.sub.2), SS(t.sub.2), S.sub.NO*(t.sub.1)), f.sub.2(S.sub.O(t.sub.2), S.sub.NH(t.sub.2), SS(t.sub.2), S.sub.NO*(t.sub.1))] is the lower-level optimization problem; {circle around (2)} set the number of particle population in the upper level optimization I.sub.1, the number of the particle population in the lower level optimization I.sub.2, a maximum number of iterations in the upper level optimization N.sub.1, and a maximum number of iterations in the lower level optimization N.sub.2, I.sub.1=50, I.sub.2=50, N.sub.1=20, N.sub.2=50; {circle around (3)} introduce a single particle swarm optimization (SPSO) algorithm to optimize the upper-level optimization problem, a position and a velocity of the ith particle are:
s.sub.i(t.sub.1)=[s.sub.i,1(t.sub.1),s.sub.i,2(t.sub.1)],(5)
v.sub.i(t.sub.1)=[v.sub.i,1(t.sub.1),v.sub.i,2(t.sub.1)],(6) s.sub.i(t.sub.1) is the position of the ith particle at time t.sub.1, s.sub.i,1(t.sub.1) is the value of S.sub.NO at time t.sub.1, s.sub.i,2(t.sub.1) is the value of MLSS at time t.sub.1, v.sub.i(t.sub.1) is the velocity of the ith particle at time t.sub.1, i is the number of particles, i=1, 2, . . . , 50, an update process of s.sub.i(t.sub.1) and v.sub.i(t.sub.1) is
v.sub.i,d(t.sub.1+1)=0.7v.sub.i,d(t.sub.1)+0.72(p.sub.i,d(t.sub.1)s.sub.i,d(t.sub.1))+0.72.sub.2(g.sub.d(t.sub.1)s.sub.i,d(t.sub.1)),(7)
s.sub.i,d(t.sub.1+1)=s.sub.i,d(t.sub.1)+v.sub.i,d(t.sub.1+1),(8) where d is the space dimension, d=1, 2, v.sub.i,d(t.sub.1) is the velocity of the ith particle in the dth dimension at time t.sub.1, p.sub.i,d(t.sub.1) is the individual optimal solution of the ith particle in the dth dimension at time t.sub.1, g.sub.d(t.sub.1) is the global optimal solutions of the ith particle at time t.sub.1; {circle around (4)} if SPSO reaches a preset maximum number of evolutions N.sub.1, stop the iterative evolution process, transfer the value of S.sub.NO* to the lower level; if SPSO does not reach the preset maximum number of evolutions N.sub.1, return to step {circle around (3)}; {circle around (5)} introduce a multiobjective particle swarm optimization (MOPSO) algorithm to optimize the lower-level optimization problem, the position of the jth particle a.sub.j(t.sub.2) and the velocity of the jth particle b.sub.j(t.sub.2) can be represented as a.sub.j(t.sub.2)=[a.sub.j,1(t.sub.2), a.sub.j,2(t.sub.2), a.sub.j,3(t.sub.2), a.sub.j,4(t.sub.2)], a.sub.i,1(t.sub.2) represents the value of S.sub.O at time t.sub.2, a.sub.i,2(t.sub.2) represents the value of S.sub.NH at time t.sub.2, a.sub.i,3(t.sub.2) represents the value of SS at time t.sub.2, a.sub.i,4(t.sub.2) represents the value of S.sub.NO* at time t.sub.2, b.sub.j(t.sub.2)=[b.sub.j,1(t.sub.2), b.sub.j,2(t.sub.2), b.sub.j,3(t.sub.2), b.sub.j,4(t.sub.2)], j is the number of particles, j=1, 2, . . . , 50; during the iterative evolution process, the obtained non-dominated solutions are conserved in the external archive Z(t.sub.2), Z(t.sub.2)=[z.sub.1(t.sub.2), z.sub.2(t.sub.2), . . . , z.sub.j(t.sub.2), . . . , z.sub.50(t.sub.2)], the update rule of the external archive is:
.sub.j(t.sub.2)=z.sub.j(t.sub.2)+0.09D(z.sub.j(t.sub.2)),(9) where z.sub.j(t.sub.2) is the jth non-dominated solution at time t.sub.2 before the archive is updated, .sub.j(t.sub.2) is the jth non-dominated solution at time t.sub.2 after the archive is updated, z.sub.j(t.sub.2)=[z.sub.j,1(t.sub.2), z.sub.j,2(t.sub.2)], .sub.j(t.sub.2)=[.sub.j,1(t.sub.2), .sub.j,2(t.sub.2)], z.sub.j,1(t.sub.2) and .sub.j,1(t.sub.2) are the values of S.sub.O before and after the archive is updated, z.sub.j,2(t.sub.2) and .sub.j,2(t.sub.2) are the values of S.sub.NO before and after the archive is updated, D is the gradient descent direction; {circle around (6)} establish a multi-input-multi-output radial basis assisted model (RBSM) based on the non-dominated solutions in Z(t.sub.2): $\begin{array}{cc}{B}_j\left(t_2\right)=\underset{k=1}{\overset{8}{.Math.}}.Math..Math.o_{j,k}\left(t_2\right){}_{j,k}\left(t_2\right),& \left(10\right)\end{array}$ where B.sub.j(t.sub.2) is an output vector of RBSM, B.sub.j(t.sub.2)=[B.sub.j,1(t.sub.2), B.sub.j,2(t.sub.2)].sup.T, B.sub.j,1(t.sub.2) is a predicted value of the aeration energy at time t.sub.2, B.sub.j,2(t.sub.2) is a predicted value of the effluent quality at time t.sub.2, o.sub.j(t.sub.2)=[o.sub.j,1(t.sub.2), o.sub.j,2(t.sub.2), . . . , o.sub.j,8(t.sub.2)].sup.T are connection weights, .sub.j(t.sub.2)=[.sub.j,1(t.sub.2), .sub.j,2(t.sub.2), . . . , .sub.j,8(t.sub.2)].sup.T is an output vector of the neurons in hidden layer, the sum of the squared errors between the output of RBSM and the actual system is expressed as
e(z.sub.n(t.sub.2))=min(B.sub.n(t.sub.2)Q(t.sub.2)).sup.T(B.sub.n(t.sub.2)Q(t.sub.2)),(11) where e(z.sub.n(t.sub.2)) is the sum of the squared errors between the outputs of the nth non-dominated solution B.sub.n(t.sub.2) and the actual system Q(t.sub.2), n[1, 50], Q(t.sub.2)=[Q.sub.1(t.sub.2), Q.sub.2(t.sub.2)] is the real outputs of the aeration energy and effluent quality in the actual system, select the solution corresponding to the minimal sum of the squared error as the global optimal solution; {circle around (7)} if MOPSO reaches a preset maximum number of evolutions N.sub.2, stop the iterative evolution process and output the optimal set-points of dissolved oxygen S.sub.O*; if MOPSO does not reach the preset maximum number of evolutions N.sub.2, return to step {circle around (5)}; (4) design a controller based on predictive control strategy: {circle around (1)} define cost functions in the predictive control strategy: $\begin{array}{cc}{J}_1\left(t\right)=\frac{1}{2}.Math.\underset{q=1}{\overset{5}{.Math.}}.Math..Math.{\left(z_1\left(t\right)-y_1\left(t\right)\right)}^2+\frac{1}{2}.Math.\underset{m=1}{\overset{4}{.Math.}}.Math..Math..Math..Math.{u\left(t\right)}^T.Math..Math..Math.u\left(t\right),.Math.J_2\left(t\right)=\frac{1}{2}.Math.\underset{q=1}{\overset{5}{.Math.}}.Math..Math.{\left(z_2\left(t\right)-y_2\left(t\right)\right)}^2+\frac{1}{2}.Math.\underset{m=1}{\overset{4}{.Math.}}.Math..Math..Math..Math.{u\left(t\right)}^T.Math..Math..Math.u\left(t\right),& \left(12\right)\end{array}$ where z.sub.1(t) and z.sub.2(t) are an optimal set-points of S.sub.O* and S.sub.NO*, y.sub.1(t) and y.sub.2(t) are predicted values of S.sub.O and S.sub.NO; {circle around (2)} update control laws based on the predictive control strategy, the updated rule is:
u(t+1)=u(t)+u(t),(13) where u(t) is the control law at time t, u(t) are control variations, whose expressions are shown as:
u.sub.1(t)=J.sub.1(t)/u.sub.1(t)+J.sub.2(t)/u.sub.1(t),
u.sub.2(t)=J.sub.1(t)/u.sub.2(t)+J.sub.2(t)/u.sub.2(t),(14) where u(t) are variations of the manipulated variables oxygen transfer coefficient K.sub.La and internal recycle Q.sub.a, u(t)=[u.sub.1(t), u.sub.2(t)]; {circle around (3)} transfer the obtained K.sub.La and Q.sub.a to an actual system of WWTP; (5) inputs of a cooperative optimal control for the wastewater treatment system are K.sub.La and Q.sub.a, outputs of the wastewater treatment system are S.sub.O and S.sub.NO, the effect of the optimal control results is reflected by the daily average of the pumping energy value, the daily average of the aeration energy value, the daily average of the effluent quality value, and the tracking control results of S.sub.O and S.sub.NO.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0041] FIG. 1 shows the control system scheme of WWTP.

[0042] FIG. 2 shows the results of PE in COCS.

[0043] FIG. 3 shows the results of AE in COCS.

[0044] FIG. 4 shows the results of EQ in COCS.

[0045] FIG. 5 shows the results of S.sub.O in COCS.

[0046] FIG. 6 shows the results of S.sub.NO in COCS.

[0047] FIG. 7 is a chart showing the contribution degree index of process variables.

DETAILED DESCRIPTION

[0048] (1) The design of the control system framework for WWTP, the steps are:

[0049] 1) the control system includes the serve, the PLC control cabinet, the blower, an ORP electrode sensor, a S.sub.O sensor, a SS sensor, a temperature sensor, a pH sensor, a S.sub.NH sensor, a TN sensor, a COD sensor, where

[0050] the serve is the operation center of control system to optimize the control procedure and realize the monitoring functions;

[0051] the PLC control cabinet is used to complete the automatic detection, monitor and control of the equipment status and operation conditions;

[0052] the blower acts as an executor of the control system to realize the aeration and ventilation of aerobic tank;

[0053] an ORP electrode sensor, located in the anaerobic tank, is configured to measure the values of ORP;

[0054] a S.sub.O sensor, located in the aerobic tank, is configured to measure the values of S.sub.O;

[0055] a SS sensor, located in the aerobic tank, is configured to measure the values of SS;

[0056] a temperature sensor, located in the effluent tank, is configured to measure the values of temperature;

[0057] a pH sensor, located in the in the effluent tank, is configured to measure the values of pH;

[0058] a S.sub.NH sensor, located in the effluent tank, is configured to measure the values of S.sub.NH;

[0059] a TN sensor, located in the effluent tank, is configured to measure the values of TN;

[0060] a COD sensor, located in the effluent tank, is configured to measure the values of COD;

[0061] 2) data collection and transmission, at first, it is necessary to achieve the data acquisition through an ORP electrode sensor, a S.sub.O sensor, a SS sensor, a temperature sensor, a pH sensor, a S.sub.NH sensor, a TN sensor, a COD sensor; then it is essential to connect PLC and the sensors by using RS 485, realize the communication, complete the upload of the collected data, and transfer the collected data to servers through coordinated communication standards;

[0062] (2) the optimal controller, which comprises the establishment of the performance indices, the development of the optimization algorithm and the design of the controller, the steps are:

[0063] 1) select the related process variables of PE: S.sub.NO, MLSS, and choose the related process variables of AE and EQ: S.sub.O, SS, S.sub.NH, S.sub.NO;

[0064] 2) formulate the two-level models based on the different time scales, the upper-level model is for PE, and the lower-level models are for AE and EQ:

F.sub.1(t.sub.1)=l.sub.1(x.sub.u(t.sub.1)),(1)

f.sub.1(t.sub.2)=l.sub.2(x.sub.l(t.sub.2),x.sub.u*(t.sub.1)),

f.sub.2(t.sub.2)=l.sub.3(x.sub.l(t.sub.2),x.sub.u*(t.sub.1)),(2)

where F.sub.1(t.sub.1) is the model for PE at time t.sub.1, l.sub.1(x.sub.u(t.sub.1)) is the mapping function of PE model, f.sub.1(t.sub.2) is the model for AE at time t.sub.2, l.sub.2(x.sub.l(t.sub.2), x*.sub.u(t.sub.1)) is the mapping function of AE model, x*.sub.u(t.sub.1) is the optimal set-points of nitrate nitrogen S.sub.NO* at time t.sub.1, f.sub.2(t.sub.2) is the model for EQ at time t.sub.2, l.sub.3(x.sub.l(t.sub.2), x*.sub.u(t.sub.1)) is the mapping function of EQ model, x.sub.u(t.sub.1)=[S.sub.NO(t.sub.1), MLSS(t.sub.1)] is the input variables vector of PE at time t.sub.1, S.sub.NO(t.sub.1) is the concentration of S.sub.NO at time t.sub.1, MLSS(t.sub.1) is the concentration of MLSS at time t.sub.1, and the initial values of the two variables are [0.85, 1.56], x.sub.l(t.sub.2)=[S.sub.O(t.sub.2), SS(t.sub.2), S.sub.NH(t.sub.2)], S.sub.O(t.sub.2) is the concentration of S.sub.O at time t.sub.2, SS(t.sub.2) is the concentration of SS at time t.sub.2, S.sub.NH(t.sub.2) is the concentration of S.sub.NH at time t.sub.2, and [S.sub.O(t.sub.2), SS(t.sub.2), S.sub.NH(t.sub.2), S.sub.NO*(t.sub.1)] is the input variables vector of AE and EQ at time t.sub.2, and the initial values are [1.9, 11.6, 3.8, 0.95];

[0065] (3) design a cooperative optimization algorithm to optimize upper-level and lower-level optimization problems to obtain the optimal set-points of the control variables, where the optimization period in the upper level is two hours, in the lower level is thirty minutes, the steps are:

[0066] {circle around (1)} formulate the upper-level and lower-level problems:

Min F.sub.1(S.sub.NO(t.sub.1),MLSS(t.sub.1)),(3)

Min[f.sub.1(S.sub.O(t.sub.2),S.sub.NH(t.sub.2),SS(t.sub.2),S.sub.NO*(t.sub.1)),

f.sub.2(S.sub.O(t.sub.2),S.sub.NH(t.sub.2),SS(t.sub.2),S.sub.NO*(t.sub.1))],(4)

where Min F.sub.1(S.sub.NO(t.sub.1), MLSS(t.sub.1)) is the upper-level optimization problem, Min [f.sub.1(S.sub.O(t.sub.2), S.sub.NH(t.sub.2), SS(t.sub.2), S.sub.NO*(t.sub.1)), f.sub.2(S.sub.O(t.sub.2), S.sub.NH(t.sub.2), SS(t.sub.2), S.sub.NO*(t.sub.1))] is the lower-level optimization problem;

[0067] {circle around (2)} set the number of the particle population in the upper level optimization I.sub.1, the number of the particle population in the lower level optimization I.sub.2, the maximum number of iterations in the upper level optimization N.sub.1, and the maximum number of iterations in the lower level optimization N.sub.2, I.sub.1=50, I.sub.2=50, N.sub.1=20, N.sub.2=50;

[0068] {circle around (3)} introduce the single particle swarm optimization (SPSO) algorithm to optimize the upper-level optimization problem, the position and the velocity of the ith particle can be shown as:

s.sub.i(t.sub.1)=[s.sub.i,1(t.sub.1),s.sub.i,2(t.sub.1)],(5)

v.sub.i(t.sub.1)=[v.sub.i,1(t.sub.1),v.sub.i,2(t.sub.1)],(6)

s.sub.i(t.sub.1) is the position of the ith particle at time t.sub.1, s.sub.i,1(t.sub.1) is the value of S.sub.NO at time t.sub.1, s.sub.i,2(t.sub.1) is the value of MLSS at time t.sub.1, v.sub.i(t.sub.1) is the velocity of the ith particle at time t.sub.1, i is the number of particles, i=1, 2, . . . , 50, an update process of s.sub.i(t.sub.1) and v.sub.i(t.sub.1) are

v.sub.i,d(t.sub.1+1)=0.7v.sub.i,d(t.sub.1)+0.72(p.sub.i,d(t.sub.1)s.sub.i,d(t.sub.1))+0.72.sub.2(g.sub.d(t.sub.1)s.sub.i,d(t.sub.1)),(7)

s.sub.i,d(t.sub.1+1)=s.sub.i,d(t.sub.1)+v.sub.i,d(t.sub.1+1),(8)

where d is the space dimension, d=1, 2, v.sub.i,d(t.sub.1) is the velocity of the ith particle in the dth dimension at time t.sub.1, p.sub.i,d(t.sub.1) is the individual optimal solution of the ith particle in the dth dimension at time t.sub.1, g.sub.d(t.sub.1) is the global optimal solutions of the ith particle at time t.sub.1;

[0069] {circle around (4)} if SPSO reaches the preset maximum number of evolutions N.sub.1, stop the iterative evolution process, transfer the value of S.sub.NO* to the lower level; if SPSO does not reach the preset maximum number of evolutions N.sub.1, return to step {circle around (3)};

[0070] {circle around (5)} introduce the multiobjective particle swarm optimization (MOPSO) algorithm to optimize the lower-level optimization problem, the position of the jth particle a.sub.j(t.sub.2) and the velocity of the jth particle b.sub.j(t.sub.2) can be represented as a.sub.j(t.sub.2)=[a.sub.j,1(t.sub.2), a.sub.j,2(t.sub.2), a.sub.j,3(t.sub.2), a.sub.j,4(t.sub.2)], a.sub.i,1(t.sub.2) represents the value of S.sub.O at time t.sub.2, a.sub.i,2(t.sub.2) represents the value of S.sub.NH at time t.sub.2, a.sub.i,3(t.sub.2) represents the value of SS at time t.sub.2, a.sub.i,4(t.sub.2) represents the value of S.sub.NO* at time t.sub.2, b.sub.j(t.sub.2)=[b.sub.j,1(t.sub.2), b.sub.j,2(t.sub.2), b.sub.j,3(t.sub.2), b.sub.j,4(t.sub.2)], j is the number of particles, j=1, 2, . . . , 50; during the iterative evolution process, the obtained non-dominated solutions are conserved in the external archive Z(t.sub.2), Z(t.sub.2)=[z.sub.1(t.sub.2), z.sub.2(t.sub.2), . . . , z.sub.j(t.sub.2), . . . , z.sub.50(t.sub.2)], the update rule of the external archive is:

.sub.j(t.sub.2)=z.sub.j(t.sub.2)+0.09D(z.sub.j(t.sub.2)),(9)

where z.sub.j(t.sub.2) is the jth non-dominated solution at time t.sub.2 before the archive is updated, .sub.j(t.sub.2) is the jth non-dominated solution at time t.sub.2 after the archive is updated, z.sub.j(t.sub.2)=[z.sub.j,1(t.sub.2), z.sub.j,2(t.sub.2)], .sub.j(t.sub.2)=[.sub.j,1(t.sub.2), .sub.j,2(t.sub.2)], z.sub.j,1(t.sub.2) and .sub.j,1(t.sub.2) are the values of S.sub.O before and after the archive is updated, z.sub.j,2(t.sub.2) and .sub.j,2(t.sub.2) are the values of S.sub.NO before and after the archive is updated, D is the gradient descent direction;

[0071] {circle around (6)} establish the multi-input-multi-output radial basis assisted model (RBSM) based on the non-dominated solutions in Z(t.sub.2):

[00003]$\begin{array}{cc}{B}_j\left(t_2\right)=\underset{k=1}{\overset{8}{.Math.}}.Math..Math.o_{j,k}\left(t_2\right){}_{j,k}\left(t_2\right),& \left(10\right)\end{array}$

where B.sub.j(t.sub.2) is the output vector of RBSM, B.sub.j(t.sub.2)=[B.sub.j,1(t.sub.2), B.sub.j,2(t.sub.2)].sup.T, B.sub.j,1(t.sub.2) is the predicted value of the aeration energy at time t.sub.2, B.sub.j,2(t.sub.2) is the predicted value of the effluent quality at time t.sub.2, o.sub.j(t.sub.2)=[o.sub.j,1(t.sub.2), o.sub.j,2(t.sub.2), . . . , o.sub.j,8(t.sub.2)].sup.T are the connection weights, .sub.j(t.sub.2)=[.sub.j,1(t.sub.2), .sub.j,2(t.sub.2), . . . , .sub.j,8(t.sub.2)].sup.T is the output vector of the neurons in hidden layer, the sum of the squared errors between the output of RBSM and the actual system is expressed as

e(z.sub.n(t.sub.2))=min(B.sub.n(t.sub.2)Q(t.sub.2)).sup.T(B.sub.n(t.sub.2)Q(t.sub.2)),(11)

where e(z.sub.n(t.sub.2)) is the sum of the squared errors between the outputs of the nth non-dominated solution B.sub.n(t.sub.2) and the actual system Q(t.sub.2), n[1, 50], Q(t.sub.2)=[Q.sub.1(t.sub.2), Q.sub.2(t.sub.2)] is the real outputs of AE and EQ in the actual system, select the solution corresponding to the minimal sum of the squared error as the global optimal solution;

[0072] {circle around (7)} if MOPSO reaches the preset maximum number of evolutions N.sub.2, stop the iterative evolution process and output the optimal set-points of dissolved oxygen S.sub.O*; if MOPSO does not reach the preset maximum number of evolutions N.sub.2, return to step {circle around (5)};

[0073] (4) design the tracking control method based on the predictive control strategy:

[0074] {circle around (1)} define the cost functions in the predictive control strategy:

[00004]$\begin{array}{cc}{J}_1\left(t\right)=\frac{1}{2}.Math.\underset{q=1}{\overset{5}{.Math.}}.Math..Math.{\left(z_1\left(t\right)-y_1\left(t\right)\right)}^2+\frac{1}{2}.Math.\underset{m=1}{\overset{4}{.Math.}}.Math..Math..Math..Math.{u\left(t\right)}^T.Math..Math..Math.u\left(t\right),.Math.J_2\left(t\right)=\frac{1}{2}.Math.\underset{q=1}{\overset{5}{.Math.}}.Math..Math.{\left(z_2\left(t\right)-y_2\left(t\right)\right)}^2+\frac{1}{2}.Math.\underset{m=1}{\overset{4}{.Math.}}.Math..Math..Math..Math.{u\left(t\right)}^T.Math..Math..Math.u\left(t\right),& \left(12\right)\end{array}$

where z.sub.1(t) and z.sub.2(t) are the optimal set-points of S.sub.O* and S.sub.NO*, y.sub.1(t) and y.sub.2(t) are the predicted values of S.sub.O and S.sub.NO;

[0075] {circle around (2)} update the control laws based on the predictive control strategy, the updated rule is:

u(t+1)=u(t)+u(t),(13)

where u(t) is the control law at time t, u(t) are the control variations, whose expressions are shown as:

[00005]$\begin{array}{cc}.Math..Math.u_1\left(t\right)=\frac{1}{2}.Math.J_1\left(t\right)/u_1\left(t\right)+\frac{1}{2}.Math.J_2\left(t\right)/u_1\left(t\right),.Math..Math..Math.u_2\left(t\right)=\frac{1}{2}.Math.J_1\left(t\right)/u_2\left(t\right)+\frac{1}{2}.Math.J_2\left(t\right)/u_2\left(t\right),& \left(14\right)\end{array}$

where u(t) are the variations of the manipulated variables oxygen transfer coefficient K.sub.La and internal recycle Q.sub.a, u(t)=[u.sub.1(t), u.sub.2(t)];

[0076] {circle around (3)} transfer the obtained K.sub.La and Q.sub.a to the actual system of WWTP;

[0077] (3) the inputs of the cooperative optimal control based wastewater treatment system is K.sub.La and Q.sub.a, the outputs of the wastewater treatment system is S.sub.O and S.sub.NO, the effect of the optimal control results is reflected by the daily average of PE value, the daily average of AE value, the daily average of EQ value, and the tracking control results of S.sub.O and S.sub.NO.

[0078] The control system scheme based on COCS is shown in FIG. 1, and the optimal control results are presented in FIGS. 2-6. FIG. 2 gives the daily average PE values, X axis shows the time, and the unit is day, Y axis is the average PE value, and the unit is /KWh. FIG. 3 gives the daily average AE values, X axis shows the time, and the unit is day, Y axis is the average AE value, and the unit is /KWh. FIG. 4 gives the daily average EQ values, X axis shows the time, and the unit is day, Y axis is the average EQ value, and the unit is /m.sup.3. FIG. 5(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. 5(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. 6(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. 6(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.