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 method of designing a cooperative optimal control system for WWTP, comprising the steps: (1) providing an optimal control system that contains a center control unit, a blower, an internal flow recycle pump, sensor unit that includes a S.sub.O sensor, two S.sub.NO sensors, a SS sensor, a S.sub.NH sensor, two MLSS sensors; the center control unit is configured to receive measured data from the sensor unit, calculate optimal set-points of control variables based on the measured data, and generate control commands to optimize operation of the WWTP; the blower and the internal flow recycle pump are applied to perform the control commands transmitted by the center control unit; the sensor unit is used to measure data of process variables, the S.sub.O sensor, located in the aerobic tank, is configured to measure values of S.sub.O; the two S.sub.NO sensors, located in the anoxic and secondary sedimentation tanks, are configured to measure values of S.sub.NO; the SS sensor, located in the aerobic tank, is configured to measure values of SS; the S.sub.NH sensor, located in the secondary sedimentation tank, is configured to measure values of S.sub.NH; the two MLSS sensors, located in the anoxic and secondary sedimentation tank, are configured to measure values of MLSS; (2) adaptive calculation of the optimal set-points of control variables, rather than the fixed values by expert experience, the steps are: 1) select related process variables of PE: S.sub.NO, MLSS, and choose related process variables of AE and EQ: S.sub.O, SS, S.sub.NH, S.sub.NO; 2) formulate two-level models based on different time scales, an upper-level model is for PE, and lower-level models are for AE and EQ: $\begin{array}{cc}{F}_1\left(t_1\right)=l_1\left(x_u\left(t_1\right)\right),& \left(1\right)\\ f_1\left(t_2\right)=l_2\left(x_l\left(t_2\right),x_u^{*}\left(t_1\right)\right),f_2\left(t_2\right)=l_3\left(x_l\left(t_2\right),x_u^{*}\left(t_1\right)\right),& \left(2\right)\end{array}$ where F.sub.1(t.sub.1) is the upper-level model for PE at time t.sub.1, l.sub.1(x.sub.u(t.sub.1)) is a mapping function of PE model, f.sub.1(t.sub.2) is the lower-level model for AE 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 AE model, x*.sub.u(t.sub.1) is optimal set-points of nitrate nitrogen S.sub.NO*, f.sub.2(t.sub.2) is the lower-level model for EQ 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 EQ model, x.sub.u(t.sub.1)=[S.sub.NO(t.sub.1), MLSS(t.sub.1)] is input variables vector of PE at time t.sub.1, S.sub.NO(t.sub.1) is concentration of S.sub.NO at time t.sub.1, MLSS(t.sub.1) is concentration of MLSS at time t.sub.1, and 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 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, [S.sub.O(t.sub.2), SS(t.sub.2), S.sub.NH(t.sub.2), S.sub.NO*(t.sub.1)] is input variables vector of AE and EQ at time t.sub.2, and initial values of [S.sub.O(t.sub.2), SS(t.sub.2), S.sub.NH(t.sub.2), S.sub.NO*(t.sub.1)] are [1.9, 11.6, 3.8, 0.95]; 3) design a cooperative optimization algorithm to optimize the upper-level and lower-level optimization problems for obtaining the optimal set-points of the control variables, where the optimization period in the upper level is 1-3 hours dedicated to the slower process dynamics, in the lower level is 10-50 minutes dedicated to the faster process dynamics, the steps are: {circle around (1)} formulate the upper-level and lower-level problems: $\begin{array}{cc}\mathrm{Min}{F}_1\left(S_{\mathrm{NO}}\left(t_1\right),\mathrm{MLSS}\left(t_1\right)\right),& \left(3\right)\\ \mathrm{Min}\left[f_1\left(S_O\left(t_2\right),S_{\mathrm{NH}}\left(t_2\right),\mathrm{SS}\left(t_2\right),S_{\mathrm{NO}}^{*}\left(t_1\right)\right),f_2\left(S_O\left(t_2\right),S_{\mathrm{NH}}\left(t_2\right),\mathrm{SS}\left(t_2\right),S_{\mathrm{NO}}^{*}\left(t_1\right)\right)\right],& \left(4\right)\end{array}$ 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 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, where I.sub.1=50, I.sub.2=50, N.sub.1=20, N.sub.2=50; {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: $\begin{array}{cc}{s}_i\left(t_1\right)=\left[s_{i,1}\left(t_1\right),s_{i,2}\left(t_1\right)\right],& \left(5\right)\\ v_i\left(t_1\right)=\left[v_{i,1}\left(t_1\right),v_{i,2}\left(t_1\right)\right],& \left(6\right)\end{array}$ 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, 1=1, 2, . . . , 50, the update process of s.sub.i(t.sub.1) and v.sub.i(t.sub.1) are $\begin{array}{cc}{v}_{i,d}\left(t_1+1\right)=0.7v_{i,d}\left(t_1\right)+0.72\left(p_{i,d}\left(t_1\right)-s_{i,d}\left(t_1\right)\right)+0.72\left(g_d\left(t_1\right)-s_{i,d}\left(t_1\right)\right),& \left(7\right)\\ s_{i,d}\left(t_1+1\right)=s_{i,d}\left(t_1\right)+v_{i,d}\left(t_1+1\right),& \left(8\right)\end{array}$ 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 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)}; {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: $\begin{array}{cc}{\stackrel{\bigvee }{z}}_j\left(t_2\right)=z_j\left(t_2\right)+0.09\nabla D\left(z_j\left(t_2\right)\right),& \left(9\right)\end{array}$ 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 the 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.}}{o}_{j,k}\left(t_2\right)\times {\theta }_{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 $\begin{array}{cc}e\left(z_n\left(t_2\right)\right)={\min \left(B_n\left(t_2\right)-Q\left(t_2\right)\right)}^T\left(B_n\left(t_2\right)-Q\left(t_2\right)\right),& \left(11\right)\end{array}$ where e(z.sub.n(t.sub.2)) is the sum of the squared errors between the outputs of 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; {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)}; then the optimal set-points of S.sub.NO*and S.sub.O* can be obtained; (3) perform the control commands: if the measured value of S.sub.NO is lower or higher than S.sub.NO*, it is noted that the internal recycle flow should be adjusted to satisfy the operating requirement by manipulating the electromagnetic valve of internal flow recycle pump; if the measured value of S.sub.O is lower or higher than S.sub.O*, it is noted that the supplied oxygen should be adjusted by manipulating the fan frequency of blower; the detailed adjusting strategy is realized by the predictive control; the steps are: {circle around (1)} define the cost functions in the predictive control strategy: $\begin{array}{cc}{J}_1\left(t\right)=\frac{1}{2}\underset{q=1}{\overset{5}{.Math.}}{\left(z_1\left(t\right)-y_1\left(t\right)\right)}^2+\frac{1}{2}\underset{m=1}{\overset{4}{.Math.}}\Delta {u\left(t\right)}^T\Delta u\left(t\right),J_2\left(t\right)=\frac{1}{2}\underset{q=1}{\overset{5}{.Math.}}{\left(z_2\left(t\right)-y_2\left(t\right)\right)}^2+\frac{1}{2}\underset{m=1}{\overset{4}{.Math.}}\Delta {u\left(t\right)}^T\Delta 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; {circle around (2)} update the control laws based on the predictive control strategy, the updated rule is: $\begin{array}{cc}{J}_1\left(t\right)=\frac{1}{2}\underset{q=1}{\overset{5}{.Math.}}{\left(z_1\left(t\right)-y_1\left(t\right)\right)}^2+\frac{1}{2}\underset{m=1}{\overset{4}{.Math.}}\Delta {u\left(t\right)}^T\Delta u\left(t\right),J_2\left(t\right)=\frac{1}{2}\underset{q=1}{\overset{5}{.Math.}}{\left(z_2\left(t\right)-y_2\left(t\right)\right)}^2+\frac{1}{2}\underset{m=1}{\overset{4}{.Math.}}\Delta {u\left(t\right)}^T\Delta u\left(t\right),& \left(12\right)\end{array}$ where u(t) is the control law at time t, Δu(t) are the control variations, whose expressions are shown as: $\begin{array}{cc}\Delta u_1\left(t\right)=\frac{1}{2}\partial J_1\left(t\right)/\partial u_1\left(t\right)+\frac{1}{2}\partial J_2\left(t\right)/\partial u_1\left(t\right),\Delta u_2\left(t\right)=\frac{1}{2}\partial J_1\left(t\right)/\partial u_2\left(t\right)+\frac{1}{2}\partial J_2\left(t\right)/\partial u_2\left(t\right),& \left(14\right)\end{array}$ where Δu(t) are the variations of the manipulated variables electromagnetic valve of internal flow recycle pump and the fan frequency of blower, Δu(t)=[Δu.sub.1(t), Δu.sub.2(t)]; the values of S.sub.O and S.sub.NO will be changed accordingly, and then transmitted to the center control unit to realize the optimal control; the effects of the proposed optimal control results are 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.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

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

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

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

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

[0034] FIG. 5 includes FIG. 5(a) and FIG. 5(b), and shows the results of S.sub.O in COCS.

[0035] FIG. 6 includes FIG. 6(a) and FIG. 6(b), and shows the results of S.sub.NO in COCS.

DETAILED DESCRIPTION

[0036] The system shown in FIG. 1 includes a sand filter for filtering out sand and other particles contained in the wastewater to be treated. Connected to the sand filter is a primary settling tank for settling filtered wastewater. The wastewater output from the primary settling tank is sent to a biochemical reaction tank to oxide and decompose the organic matters. After reaction in a biochemical reaction tank, the wastewater is introduced into a secondary sedimentation tank, after sedimentation treatment, the wastewater is discharged to complete the wastewater treatment.

[0037] 1. The design of the cooperative optimal control system for WWTP, the steps are:

[0038] (1) The designed optimal control system contains the center control unit, blower, internal flow recycle pump, sensor unit that includes a S.sub.O sensor, two S.sub.NO sensors, a SS sensor, a S.sub.NH sensor, two MLSS sensors.

[0039] The center control unit is used to receive the measured data from the installed sensors, calculate the optimal set-points of the control variables based on the measured data, and generate the control commands to guarantee the operation of WWTP. The blower and internal flow recycle pump are applied to perform the control commands transmitted by the center control unit. The sensor unit is used to measure the data of process variables, a S.sub.O sensor, located in the aerobic tank, is configured to measure the values of S.sub.O; two S.sub.NO sensors, located in the anoxic and secondary sedimentation tanks, are configured to measure the values of S.sub.NO; a SS sensor, located in the aerobic tank, is configured to measure the values of SS; a S.sub.NH sensor, located in the secondary sedimentation tank, is configured to measure the values of S.sub.NH; two MLSS sensors, located in the anoxic and secondary sedimentation tank, are configured to measure the values of MLSS.

[0040] (2) The adaptive calculation of the optimal set-points of control variables, rather than the fixed values by expert experience, the steps are:

[0041] 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;

[0042] 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:

[00011]$\begin{array}{cc}{F}_1\left(t_1\right)=l_1\left(x_u\left(t_1\right)\right),& \left(1\right)\\ f_1\left(t_2\right)=l_2\left(x_l\left(t_2\right),x_u^{*}\left(t_1\right)\right),f_2\left(t_2\right)=l_3\left(x_l\left(t_2\right),x_u^{*}\left(t_1\right)\right),& \left(2\right)\end{array}$

where F.sub.1(t.sub.1) is the upper-level model for PE at time t.sub.1, l.sub.1(x.sub.u(t.sub.1)) is a mapping function of PE model, f.sub.1(t.sub.2) is the lower-level model for AE 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 AE model, x*.sub.u(t.sub.1) is optimal set-points of nitrate nitrogen S.sub.NO*, f.sub.2(t.sub.2) is the lower-level model for EQ 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 EQ model, x.sub.u(t.sub.1)=[S.sub.NO(t.sub.1), MLSS(t.sub.1)] is input variables vector of PE at time t.sub.1, S.sub.NO(t.sub.1) is concentration of S.sub.NO at time t.sub.1, MLSS(t.sub.1) is concentration of MLSS at time t.sub.1, and 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 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, [S.sub.O(t.sub.2), SS(t.sub.2), S.sub.NH(t.sub.2), S.sub.NO*(t.sub.1)] is input variables vector of AE and EQ at time t.sub.2, and initial values of [S.sub.O(t.sub.2), SS(t.sub.2), S.sub.NH(t.sub.2), S.sub.NO*(t.sub.1)] are [1.9, 11.6, 3.8, 0.95];

[0043] 3) design a cooperative optimization algorithm to optimize the upper-level and lower-level optimization problems for obtaining the optimal set-points of the control variables, where the optimization period in the upper level is 1-3 hours, such as 2 hours, dedicated to the slower process dynamics, in the lower level is 10-50 minutes, such as 30 minutes, dedicated to the faster process dynamics, the steps are:

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

[00012]$\begin{array}{cc}\mathrm{Min}{F}_1\left(S_{\mathrm{NO}}\left(t_1\right),\mathrm{MLSS}\left(t_1\right)\right),& \left(3\right)\\ \mathrm{Min}\left[f_1\left(S_O\left(t_2\right),S_{\mathrm{NH}}\left(t_2\right),\mathrm{SS}\left(t_2\right),S_{\mathrm{NO}}^{*}\left(t_1\right)\right),f_2\left(S_O\left(t_2\right),S_{\mathrm{NH}}\left(t_2\right),\mathrm{SS}\left(t_2\right),S_{\mathrm{NO}}^{*}\left(t_1\right)\right)\right],& \left(4\right)\end{array}$

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;

[0045] {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, where I.sub.1=50, I.sub.2=50, N.sub.1=20, N.sub.2=50;

[0046] {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:

[00013]$\begin{array}{cc}{s}_i\left(t_1\right)=\left[s_{i,1}\left(t_1\right),s_{i,2}\left(t_1\right)\right],& \left(5\right)\\ v_i\left(t_1\right)=\left[v_{i,1}\left(t_1\right),v_{i,2}\left(t_1\right)\right],& \left(6\right)\end{array}$

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, the update process of s.sub.i(t.sub.1) and v.sub.i(t.sub.1) are

[00014]$\begin{array}{cc}{v}_{i,d}\left(t_1+1\right)=0.7v_{i,d}\left(t_1\right)+0.72\left(p_{i,d}\left(t_1\right)-s_{i,d}\left(t_1\right)\right)+0.72\left(g_d\left(t_1\right)-s_{i,d}\left(t_1\right)\right),& \left(7\right)\\ s_{i,d}\left(t_1+1\right)=s_{i,d}\left(t_1\right)+v_{i,d}\left(t_1+1\right),& \left(8\right)\end{array}$

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;

[0047] {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)};

[0048] {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:

[00015]$\begin{array}{cc}{\stackrel{\bigvee }{z}}_j\left(t_2\right)=z_j\left(t_2\right)+0.09\nabla D\left(z_j\left(t_2\right)\right),& \left(9\right)\end{array}$

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, which are derived by MOPSO algorithm, 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, which are derived by MOPSO algorithm, VD is the gradient descent direction;

[0049] {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):

[00016]$\begin{array}{cc}{B}_j\left(t_2\right)=\underset{k=1}{\overset{8}{.Math.}}{o}_{j,k}\left(t_2\right)\times {\theta }_{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

[00017]$\begin{array}{cc}e\left(z_n\left(t_2\right)\right)={\min \left(B_n\left(t_2\right)-Q\left(t_2\right)\right)}^T\left(B_n\left(t_2\right)-Q\left(t_2\right)\right),& \left(11\right)\end{array}$

where e(z.sub.n(t.sub.2)) is the sum of the squared errors between the outputs of 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 that obtained by electricity meter and the sensors in real time, select the solution corresponding to the minimal sum of the squared error as the global optimal solution;

[0050] {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)}; then the optimal set-points of S.sub.NO*and S.sub.O* can be obtained.

[0051] (3) Perform the control commands. If the measured value of S.sub.NO is lower or higher than S.sub.NO*, it is noted that the internal recycle flow in anoxic tank should be adjusted to satisfy the operating requirement by manipulating the electromagnetic valve of the internal flow recycle pump. Specifically, the signal or data of the measured value of S.sub.NO can be transmitted to the center control unit in real time, and the center control unit compares the measured value of S.sub.NO with S.sub.NO*, and then sends a command to the internal flow recycle pump to control its electromagnetic valve in real time. If the measured value of S.sub.O is lower or higher than S.sub.O*, it is noted that the supplied oxygen should be adjusted by manipulating the fan frequency of the blower. Specifically, the signal or data of the measured value of S.sub.O can be transmitted to the center control unit in real time, and the center control unit compares the measured value of S.sub.O with S.sub.O*, and then sends a command to the blower to control its fan frequency in real time. The detailed adjusting strategy is realized by the predictive control. The steps are:

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

[00018]$\begin{array}{cc}{J}_1\left(t\right)=\frac{1}{2}\underset{q=1}{\overset{5}{.Math.}}{\left(z_1\left(t\right)-y_1\left(t\right)\right)}^2+\frac{1}{2}\underset{m=1}{\overset{4}{.Math.}}\Delta {u\left(t\right)}^T\Delta u\left(t\right),J_2\left(t\right)=\frac{1}{2}\underset{q=1}{\overset{5}{.Math.}}{\left(z_2\left(t\right)-y_2\left(t\right)\right)}^2+\frac{1}{2}\underset{m=1}{\overset{4}{.Math.}}\Delta {u\left(t\right)}^T\Delta 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;

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

[00019]$\begin{array}{cc}u\left(t+1\right)=u\left(t\right)+\Delta u\left(t\right),& \left(13\right)\end{array}$

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

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

where Δu(t) are the variations of the manipulated variables electromagnetic valve of internal flow recycle pump and the fan frequency of blower, Δu(t)=[Δu.sub.1(t), Δu.sub.2(t)]. The values of S.sub.O and S.sub.NO will be changed accordingly, and then transmitted to the center control unit to realize the optimal control. The effects of the proposed optimal control results are 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.

[0054] 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.