Data-knowledge driven optimal control method for municipal wastewater treatment process

Abstract

A data-knowledge driven multi-objective optimal control method for municipal wastewater treatment process belongs to the field of wastewater treatment. To balance the energy consumption and effluent quality, a data driven multi-objective optimization model, including energy consumption model and effluent quality model are established to obtain the nonlinear relationship along energy consumption, effluent quality and manipulated variables. Meanwhile, a multi-objective particle swarm optimization algorithm, based on evolutionary knowledge, is proposed to optimize the set-points of nitrate nitrogen and dissolved oxygen. Moreover, the proportional integral differential (PID) controller is designed to track the set-points. Then the effluent quality can be improved and the energy consumption can be reduced.

Claims

1. A data-knowledge driven optimal control method of municipal wastewater treatment process, comprising obtaining optimal set-points of manipulated variables and tracking the manipulated variables to improve effluent quality and reduce energy consumption, the method comprising the following technical scheme and implementation steps: (1) establish data-driven multi-objective optimization model: I. taking energy consumption and effluent quality as objectives, a multi-objective optimization model is established for municipal wastewater treatment process.
min F(t)=[ƒ.sub.1(t),ƒ.sub.2(t)]  (1)  where F(t) is the multi-objective optimization model of municipal wastewater treatment process at time t, ƒ.sub.1(t) is energy consumption at time t, ƒ.sub.2(t) is effluent quality at time t; II. data-driven energy consumption and effluent quality models are established as $\begin{array}{cc}{f}_1\left(t\right)=W_{10}\left(t\right)+\underset{i=1}{\overset{I_1}{.Math.}}{W}_{1i}\left(t\right)B_{1i}\left(t\right)& \left(2\right)\\ f_2\left(t\right)=W_{20}\left(t\right)+\underset{i=1}{\overset{I_2}{.Math.}}{W}_{2i}\left(t\right)B_{2i}\left(t\right)& \left(3\right)\end{array}$  where I.sub.1 is a number of radial basis kernel functions of the energy consumption model, I.sub.1∈[3, 30], I.sub.2 is a number of radial basis kernel functions of the effluent quality model, I.sub.2∈[3, 30], W.sub.10(t) is an output offset of the energy consumption model, W.sub.20(t) is an output offset of the effluent quality model, W.sub.1i(t) is a weight of ith radial basis kernel function in the energy consumption model, W.sub.2i(t) is a weight of ith radial basis kernel function in the effluent quality model, B.sub.1i(t) is ith radial basis kernel function related to the energy consumption model, B.sub.2i(t) is ith radial basis kernel function related to the effluent quality model. $\begin{array}{cc}{B}_{1i}\left(t\right)=e^{-{.Math.s\left(t\right)-c_{1i}\left(t\right).Math.}^2/2{{\sigma }_{1i}\left(t\right)}^2}& \left(4\right)\\ B_{2i}\left(t\right)=e^{-{.Math.s\left(t\right)-c_{2i}\left(t\right).Math.}^2/2{{\sigma }_{2i}\left(t\right)}^2}& \left(5\right)\end{array}$  where s(t)=[S.sub.NO(t), S.sub.O(t), MLSS(t), S.sub.NH(t)] is an input vector, S.sub.NO(t) is a concentration of nitrate nitrogen in anaerobic final stage at time t, S.sub.NO(t)∈[0.2 mg/L, 2 mg/L], S.sub.O(t) is a concentration of dissolved oxygen in aerobic end stage at time t, S.sub.O(t)∈[0.4 mg/L, 3 mg/L], MLSS(t) is an effluent concentration of mixed liquor suspended solids at time t, MLSS(t)∈[0 mg/L, 100 mg/L], S.sub.NH(t) is an effluent concentration of ammonia nitrogen at time t, S.sub.NH(t)∈[0 mg/L, 4 mg/L], c.sub.1i(t) is a center of ith radial basis function in the energy consumption model, all variables of c.sub.1i(t) are limited in [−1, 1], c.sub.2i(t) is a center of ith radial basis function in the effluent quality model, all variables of c.sub.2i(t) are limited in [−1, 1], σ.sub.1i(t) is a width of ith radial basis function in the energy consumption model, σ.sub.1i(t)∈[0, 3], σ.sub.2i(t) is a width of ith radial basis function in the effluent quality model, σ.sub.2i(t)∈[0, 3]; (2) design multi-objective particle swarm optimization based on evolutionary knowledge: 1) controllable variables S.sub.NO and S.sub.O of municipal wastewater treatment process are used as position variables of multi-objective particle swarm optimization; population size of multi-objective particle swarm optimization is set to N, N∈[10, 100]; maximum iteration time of multi-objective particle swarm optimization is set to K, K∈[50, 200]; iteration time of population is set to k, k∈[1, K]; a number of iterations of particle information is set to k.sub.0, k.sub.0∈[2, 10]; 2) initialize population: a population with N particles is randomly generated; objective values are obtained by formula (1); a personal best position is
p.sub.n(1)=x.sub.n(1)  (6)  where p.sub.n(1) is the personal best position of nth particle in a first iteration, x.sub.n(1)=[x.sub.n,1(1), x.sub.n,2(1)] is a position of the nth particle in the first iteration, x.sub.n,1(1) is a first dimensional position of the nth particle in the first iteration, x.sub.n,1(1)∈[0.2 mg/L, 2 mg/L], x.sub.n,2(1) is a second dimensional position of the nth particle in the first iteration, x.sub.n,2(1)∈[0.4 mg/L, 3 mg/L]; establish archive A(1): an archive is obtained by comparing objectives between particles; when both objectives of a particle are less than or equal to corresponding objectives of other particles, and at least one objective is smaller than the corresponding objective of other particles, then this particle is called a non-dominated solution; by comparing the objectives of particle, the non-dominated solutions are stored in the archive; a diversity distribution is calculated by $\begin{array}{cc}D{S}_n\left(1\right)=\underset{m=1}{\overset{M}{.Math.}}\underset{j=1}{\overset{N}{.Math.}}.Math.(f_{n,m}\left(1\right)-f_{j,m}\left(1\right).Math.\text{/}N& \left(7\right)\end{array}$  where DS.sub.n(1) is a diversity distribution of the nth particle in the first iteration, ƒ.sub.n,m(1) is mth objective value of the nth particle in the first iteration, |⋅| represents absolute value; 3) evolutionary process of population I. enter next iteration by increasing the number of iterations by 1; a convergence distribution and a diversity distribution of each particle are recorded in the evolutionary process of population: $\begin{array}{cc}C{S}_n\left(k\right)=\{\begin{array}{cc}\underset{m=1}{\overset{M}{.Math.}}\left(f_{n,m}\left(k-1\right)-f_{n,m}\left(k\right)\right),& \mathrm{if}{x}_n\left(k\right)<x_n\left(k-1\right)\\ 0,& \mathrm{otherwhise}\end{array}& \left(8\right)\\ {\mathrm{DS}}_n\left(k\right)=\underset{m=1}{\overset{M}{.Math.}}\underset{j=1}{\overset{N}{.Math.}}.Math.(f_{n,m}\left(k\right)-f_{j,m}\left(k\right).Math.\text{/}N& \left(9\right)\end{array}$  where CS.sub.n(k) is a convergence distribution of nth particle in kth iteration, ƒ.sub.n,m(k) is mth objective value of the nth particle in the kth iteration, m∈[1, M], M=2, x.sub.n(k) is a position vector of the nth particle, DS.sub.n(k) is a diversity distribution of the nth particle in the kth iteration, |⋅| represents absolute value; II. convergence and diversity indexes of individual and population are established by using distribution knowledge, in which the distribution knowledge includes historical distributions of particles. $\begin{array}{cc}{\mathrm{IC}}_n\left(k\right)=\underset{u=k-k_0}{\overset{k}{.Math.}}{e}^{-\frac{C{S}_n\left(k\right)}{k-u+1}}& \left(10\right)\\ \mathrm{PC}\left(k\right)=\underset{n=1}{\overset{N}{.Math.}}{\mathrm{IC}}_n\left(k\right)& \left(11\right)\\ {ID}_n\left(k\right)=\underset{u=k-k_0}{\overset{k}{.Math.}}{e}^{-\frac{D{S}_n\left(k\right)}{k-u+1}}& \left(12\right)\\ \mathrm{PD}\left(k\right)=\underset{n=1}{\overset{N}{.Math.}}{\mathrm{ID}}_n\left(k\right)& \left(13\right)\end{array}$  where IC.sub.n(k) is individual convergence of the nth particle in the kth iteration, PC(k) is population convergence in the kth iteration, ID.sub.n(k) is individual diversity of the nth particle in the kth iteration, PD(k) is population diversity in the kth iteration, u∈[k−k.sub.0, k] is iteration times; III. select evolutionary strategy of population: Case 1: when PC(k)>PC(k−1) and PD(k)>PD(k−1), velocity and position of particle are updated by
v.sub.n,d(k+1)=ωv.sub.n,d(k)+c.sub.1r.sub.1(p.sub.n,d(k)−x.sub.n,d(k))+c.sub.2r.sub.2(g.sub.d(k)−x.sub.i,d(k))  (14)
x.sub.n,d(k+1)=x.sub.n,d(k)+v.sub.n,d(k+1)  (15)  where ω is inertia weight selected in [0.5, 0.9] randomly, v.sub.n,d(k) is d-dimensional velocity of the nth particle in the kth iteration, x.sub.n,d(k) is d-dimensional position of the nth particle in the kth iteration, p.sub.n,d(k) is d-dimensional personal best position of the nth particle in the kth iteration, g.sub.d(k) is d-dimensional position of the population in the kth iteration, r.sub.1 and r.sub.2 are a random value distributed in [0, 1], c.sub.1 is an acceleration factor of personal best solution, selected in [1.5, 2.5] randomly, c.sub.2 is an acceleration factor of global best solution, selected in [1.5, 2.5] randomly; Case 2: when PC(k)<PC(k−1) and PD(k)>PD(k−1), the velocity and position of particle are updated by
v.sub.n,d(k+1)=ωv.sub.n,d(k)+c.sub.1r.sub.1(p.sub.n,d(k)−x.sub.n,d(k))+c.sub.2r.sub.2(g.sub.d(k)−x.sub.i,d(k))+c.sub.3r.sub.3C.sub.d(k)  (16)
x.sub.n,d(k+1)=x.sub.n,d(k)+v.sub.n,d(k+1)  (17)  where r.sub.3 is a random value distributed in [0, 1], c.sub.3 is an acceleration factor related to convergence direction, selected in [0.3, 0.5] randomly, C.sub.d(k) is d-dimensional flight direction of particles with maximum convergence in the population; Case 3: when PC(k)>PC(k−1) and PD(k)<PD(k−1), the velocity and position of particle are updated by
v.sub.n,d(k+1)=ωv.sub.n,d(k)+c.sub.1r.sub.1(p.sub.n,d(k)−x.sub.n,d(k))+c.sub.2r.sub.2(g.sub.d(k)−x.sub.i,d(k))+c.sub.4r.sub.4D.sub.d(k)  (18)
x.sub.n,d(k+1)=x.sub.n,d(k)+v.sub.n,d(k+1)  (19)  where r.sub.4 is a random value distributed in [0, 1], c.sub.4 is an acceleration factor related to diversity direction, selected in [0.3, 0.5] randomly, D.sub.d(k) is d-dimensional flight direction of particles with maximum diversity in the population; Case 4: when PC(k)<PC(k−1) and PD(k)<PD(k−1), the velocity and position of particle are updated by
v.sub.n,d(k+1)=ωv.sub.n,d(k)+c.sub.1r.sub.1(p.sub.n,d(k)−x.sub.n,d(k))+c.sub.2r.sub.2(g.sub.d(k)−x.sub.i,d(k))+½(c.sub.3r.sub.3C.sub.d(k)+c.sub.4r.sub.4D.sub.d(k))   (20)
x.sub.n,d(k+1)=x.sub.n,d(k)+v.sub.n,d(k+1)  (21) Case 5: when PC(k)=PC(k−1) or PD(k)=PD(k−1), the velocity and position of particle are updated by $\begin{array}{cc}{v}_{n,d}\left(k+1\right)=\omega v_{n,d}\left(k\right)+c_1{r}_1\left(p_{n,d}\left(k\right)-x_{n,d}\left(k\right)\right)+c_2{r}_2\left(g_d\left(k\right)-x_{i,d}\left(k\right)\right)& \left(22\right)\\ x_{n,d}\left(k+1\right)=\{\begin{array}{cc}{x}_{d,\min }+\left(x_{d,\max }-x_{d,\min }\right)\times U\left(0,1\right),& r_5\le p_b\\ x_{n,d}\left(k\right),& r_5>p_b\end{array}& \left(23\right)\end{array}$  where U(0, 1) is a random value between 0 and 1 which obeys uniform distribution, x.sub.d,min is a minimum boundary value of d-dimensional position, x.sub.d,max is a maximum boundary value of d-dimensional position; when d=1, x.sub.1,min=0.2 mg/L, x.sub.1,max=2 mg/L; when d=2, x.sub.2,min=0.4 mg/L, x.sub.2,max=3 mg/L; r.sub.5 is a random value distributed in [0, 1], p.sub.b is mutation probability, which is described as $\begin{array}{cc}{p}_b=0.5-0.5\times \frac{k}{K}& \left(24\right)\end{array}$ IV. population in the kth iteration is combined with archive A(k−1) to obtain J(k), and then the non-dominated solutions are selected from J(k) to establish A(k); V. if k is greater than or equal to K, go to step VI; if k is less than K, go to step I; VI. in the archive A(K), a non-dominated solution is randomly selected as an optimal set-point a*(t)=a.sub.h(K), a.sub.h(K)=[S.sub.NO*(K), S.sub.O*(K)], where S.sub.NO*(K) is an optimal set-point of S.sub.NO, S.sub.O*(K) is an optimal set-point of S.sub.O; then, the optimal set-points are saved; (3) track the optimal set-points S.sub.NO*(K) and S.sub.O*(K): PID controller is designed to track the optimal set-points S.sub.NO*(K) and S.sub.O*(K); the expression of PID controller is $\begin{array}{cc}\Delta z\left(t\right)=K_p\left[e\left(t\right)+H_l{\int }_0^{t}e\left(t\right)dt+H_d\frac{de\left(t\right)}{dt}\right]& \left(25\right)\\ K_p=\left[\begin{array}{cc}10000& 0\\ 0& 20\end{array}\right]& \left(26\right)\\ H_l=\left[\begin{array}{cc}3000& 0\\ 0& 5\end{array}\right]& \left(27\right)\\ H_d=\left[\begin{array}{cc}100& 0\\ 0& 1\end{array}\right]& \left(28\right)\end{array}$  where Δz(t)=[ΔQ.sub.a(t), ΔK.sub.La5(t)].sup.T is a manipulated variable matrix, ΔQ.sub.a(t) is a change of internal circulation flow in wastewater treatment process, ΔK.sub.La5(t) is a change of oxygen transfer coefficient in a fifth zone of biochemical reactor, K.sub.p is a proportional coefficient matrix, H.sub.l is an integral coefficient matrix, and H.sub.d is a differential coefficient matrix; e(t)=y*(t).sup.T−y(t).sup.T is a control error, y*(t)=[S.sub.NO*(t), S.sub.O*(t)] is an optimal set-point matrix at time t, y(t)=[S.sub.NO(t), S.sub.O(t)] is an actual output matrix; (4) inputs of the data-knowledge driven optimal control method of municipal wastewater treatment process are the change of internal circulation flow ΔQ.sub.a(t) and the change of oxygen transfer coefficient in the fifth zone of biochemical reactor ΔK.sub.La5(t); the optimal set-points of S.sub.NO and S.sub.O in municipal wastewater treatment process are tracked and controlled.

Description

DESCRIPTION OF DRAWINGS

[0033] FIG. 1 shows the framework of data-knowledge driven optimal control method.

[0034] FIG. 2 shows the tracking result of nitrate nitrogen for the optimal control method.

[0035] FIG. 3 shows the tracking error of nitrate nitrogen for the optimal control method.

[0036] FIG. 4 shows the tracking result of dissolved oxygen for the optimal control method.

[0037] FIG. 5 shows the tracking error of dissolved oxygen for the optimal control method.

DETAILED DESCRIPTION OF THE INVENTION

[0038] A data-knowledge driven optimal control method is designed for municipal wastewater treatment process in this patent. Its characteristic lies in obtaining the optimal set-points of manipulated variables and tracking these variables to improve effluent quality and reduce energy consumption. This patent adopts the following technical scheme and implementation steps: [0039] (1) Establish data-driven multi-objective optimization model:

[0040] I. Taking energy consumption and effluent quality as objectives, a multi-objective optimization model is established for municipal wastewater treatment process.

min F(t)=[ƒ.sub.1(t),ƒ.sub.2(t)]  (1)

where F(t) is the multi-objective optimization model of municipal wastewater treatment process at time t, ƒ.sub.1(t) is the energy consumption at time t, ƒ.sub.2(t) is the effluent quality at time t;

[0041] II. The data-driven energy consumption and effluent quality models are established as

[00009]$\begin{array}{cc}{f}_1\left(t\right)=W_{10}\left(t\right)+\underset{i=1}{\overset{I_1}{.Math.}}{W}_{1i}\left(t\right)B_{1i}\left(t\right)& \left(2\right)\\ f_2\left(t\right)=W_{20}\left(t\right)+\underset{i=1}{\overset{I_2}{.Math.}}{W}_{2i}\left(t\right)B_{2i}\left(t\right)& \left(3\right)\end{array}$

where I.sub.1=10 is the number of radial basis kernel functions of energy consumption model, I.sub.2=10 is the number of radial basis kernel functions of effluent quality model, W.sub.10(t)=−1.20 is the output offset of energy consumption model, W.sub.10(0)=−1.20, W.sub.20(t) is the output offset of effluent quality model, W.sub.20(0)=0.34, W.sub.1i(t) is the weight of the ith radial basis kernel function in energy consumption model, W.sub.1i(0)=−0.78, W.sub.2i(t) is the weight of the ith radial basis kernel function in effluent quality model, W.sub.2i(0)=1.62, B.sub.1i(t) is the ith radial basis kernel function related to energy consumption model, B.sub.2i(t) is the ith radial basis kernel function related to effluent quality model.

[00010]$\begin{array}{cc}{B}_{1i}\left(t\right)=e^{-{.Math.s\left(t\right)-c_{1i}\left(t\right).Math.}^2/2{{\sigma }_{1i}\left(t\right)}^2}& \left(4\right)\\ B_{2i}\left(t\right)=e^{-{.Math.s\left(t\right)-c_{2i}\left(t\right).Math.}^2/2{{\sigma }_{2i}\left(t\right)}^2}& \left(5\right)\end{array}$

where s(t)=[S.sub.NO(t), S.sub.O(t), MLSS(t), S.sub.NH(t)] is the input vector, s(0)=[1, 1.5, 15, 2.3], S.sub.NO(t) is the concentration of nitrate nitrogen in anaerobic final stage at time t, S.sub.NO(t)∈[0.2 mg/L, 2 mg/L], S.sub.O(t) is the concentration of dissolved oxygen in aerobic end stage at time t, S.sub.O(t)∈[0.4 mg/L, 3 mg/L], MLSS(t) is the effluent concentration of mixed liquor suspended solids at time t, MLSS(t)∈[0 mg/L, 100 mg/L], S.sub.NH(t) is the effluent concentration of ammonia nitrogen at time t, S.sub.NH(t)∈[0 mg/L, 4 mg/L], c.sub.1i(t) is the center of the ith radial basis function in energy consumption model, c.sub.1i(0)=[0.76, 0.45, 0.21, −0.33], c.sub.2i(t) is the center of the ith radial basis function in effluent quality model, c.sub.2i(0)=[0.82, 0.67, −0.29, 0.85], σ.sub.1i(t) is the width of the ith radial basis function in the energy consumption model, σ.sub.2i(0)=0.62, σ.sub.2i(t) is the width of the ith radial basis function in the effluent quality model, σ.sub.2i(0)=1.72;

[0042] (2) Design multi-objective particle swarm optimization based on evolutionary knowledge:

[0043] 1) The controllable variables S.sub.NO and S.sub.O of municipal wastewater treatment process are used as the position variables of multi-objective particle swarm optimization. The population size of multi-objective particle swarm optimization is set to N, N=20. The maximum iteration time of multi-objective particle swarm optimization is set to K, K=100. The iteration time of population is set to k, k∈[1, K]. The number of iterations of particle information is set to k.sub.0=4;

[0044] 2) Initialize the population: the population with N particles is randomly generated. The objective values are obtained by formula (1). The personal best position is

p.sub.n(1)=x.sub.n(1)  (6)

where p.sub.n(1) is the personal best position of the nth particle in the first iteration, x.sub.n(1)=[x.sub.n,1(1), x.sub.n,2(1)] is the position of the nth particle in the first iteration, x.sub.n,1(1) is the first dimensional position of the nth particle in the first iteration, x.sub.n,1(1)∈[0.2 mg/L, 2 mg/L], x.sub.n,2(1) is the second dimensional position of the nth particle in the first iteration, x.sub.n,2(1)∈[0.4 mg/L, 3 mg/L];

[0045] Establish the archive A(1): the archive is obtained by comparing the objectives between particles. When both objectives of a particle are less than or equal to the corresponding objectives of other particles, and at least one objective is smaller than the corresponding objective of other particles, then this particle is called the non-dominated solution. By comparing the objectives of particle, the non-dominated solutions are stored in the archive;

[0046] The diversity distribution is calculated by

[00011]$\begin{array}{cc}D{S}_n\left(1\right)=\underset{m=1}{\overset{M}{.Math.}}\underset{j=1}{\overset{N}{.Math.}}.Math.(f_{n,m}\left(1\right)-f_{j,m}\left(1\right).Math.\text{/}N& \left(7\right)\end{array}$

where DS.sub.n(1) is the diversity distribution of the nth particle in the first iteration, ƒ.sub.n,m(1) is the mth objective value of the nth particle in the first iteration, |⋅| represents absolute value;

[0047] 3) The evolutionary process of population

[0048] I. Enter the next iteration, that is, increase the number of iterations by 1. The convergence distribution and diversity distribution of each particle are recorded in the evolutionary process:

[00012]$\begin{array}{cc}C{S}_n\left(k\right)=\{\begin{array}{cc}\underset{m=1}{\overset{M}{.Math.}}\left(f_{n,m}\left(k-1\right)-f_{n,m}\left(k\right)\right),& \mathrm{if}{x}_n\left(k\right)<x_n\left(k-1\right)\\ 0,& \mathrm{otherwhise}\end{array}& \left(8\right)\\ {\mathrm{DS}}_n\left(k\right)=\underset{m=1}{\overset{M}{.Math.}}\underset{j=1}{\overset{N}{.Math.}}.Math.(f_{n,m}\left(k\right)-f_{j,m}\left(k\right).Math.\text{/}N& \left(9\right)\end{array}$

where CS.sub.n(k) is the convergence distribution of the nth particle in the kth iteration, ƒ.sub.n,m(k) is the mth objective value of the nth particle in the kth iteration, m∈[1, M], M=2, x.sub.n(k) is the position vector of the nth particle, DS.sub.n(k) is the diversity distribution of the nth particle in the kth iteration, |⋅| represents absolute value;

[0049] II. The convergence and diversity indexes of individual and population are established by using distribution knowledge, in which the distribution knowledge consists of historical distributions of particles.

[00013]$\begin{array}{cc}{\mathrm{IC}}_n\left(k\right)=\underset{u=k-k_0}{\overset{k}{.Math.}}{e}^{-\frac{C{S}_n\left(k\right)}{k-u+1}}& \left(10\right)\\ \mathrm{PC}\left(k\right)=\underset{n=1}{\overset{N}{.Math.}}{\mathrm{IC}}_n\left(k\right)& \left(11\right)\\ {ID}_n\left(k\right)=\underset{u=k-k_0}{\overset{k}{.Math.}}{e}^{-\frac{D{S}_n\left(k\right)}{k-u+1}}& \left(12\right)\\ \mathrm{PD}\left(k\right)=\underset{n=1}{\overset{N}{.Math.}}{\mathrm{ID}}_n\left(k\right)& \left(13\right)\end{array}$

where IC.sub.n(k) is the individual convergence of the nth particle in the kth iteration, PC(k) is the population convergence in the kth iteration, ID.sub.n(k) is the individual diversity of the nth particle in the kth iteration, PD(k) is the population diversity in the kth iteration, u∈[k−k.sub.0, k] is the iteration times;

[0050] III. Select the evolutionary strategy of population:

[0051] Case 1: When PC(k)>PC(k−1) and PD(k)>PD(k−1), the velocity and position of particle are updated by

v.sub.n,d(k+1)=ωv.sub.n,d(k)+c.sub.1r.sub.1(p.sub.n,d(k)−x.sub.n,d(k))+c.sub.2r.sub.2(g.sub.d(k)−(k))  (14)

x.sub.n,d(k+1)=x.sub.n,d(k)+v.sub.n,d(k+1)  (15)

where ω is the inertia weight selected in [0.5, 0.9] randomly, v.sub.n,d(k) is the d-dimensional velocity of the nth particle in the kth iteration, x.sub.n,d(k) is the d-dimensional position of the nth particle in the kth iteration, p.sub.n,d(k) is the d-dimensional personal best position of the nth particle in the kth iteration, g.sub.d(k) is the d-dimensional position of the population in the kth iteration, r.sub.1 and r.sub.2 are the random value distributed in [0, 1], c.sub.1 is the acceleration factor of personal best solution, selected in [1.5, 2.5] randomly, c.sub.2 is the acceleration factor of global best solution, selected in [1.5, 2.5] randomly.

[0052] Case 2: When PC(k)<PC(k−1) and PD(k)>PD(k−1), the velocity and position of particle are updated by

v.sub.n,d(k+1)=ωv.sub.n,d(k)+c.sub.1r.sub.1(p.sub.n,d(k)−x.sub.n,d(k))+c.sub.2r.sub.2(g.sub.d(k)−x.sub.i,d(k))+c.sub.3r.sub.3C.sub.d(k)  (16)

x.sub.n,d(k+1)=x.sub.n,d(k)+v.sub.n,d(k+1)  (17)

where r.sub.3 is the random value distributed in [0, 1], c.sub.3 is the acceleration factor related to convergence direction, selected in [0.3, 0.5] randomly, C.sub.d(k) is the d-dimensional flight direction of particles with maximum convergence in the population.

[0053] Case 3: When PC(k)>PC(k−1) and PD(k)<PD(k−1), the velocity and position of particle are updated by

v.sub.n,d(k+1)=ωv.sub.n,d(k)+c.sub.1r.sub.1(p.sub.n,d(k)−x.sub.n,d(k))+c.sub.2r.sub.2(g.sub.d(k)−x.sub.i,d(k))+c.sub.4r.sub.4D.sub.d(k)  (18)

x.sub.n,d(k+1)=x.sub.n,d(k)+v.sub.n,d(k+1)  (19)

where r.sub.4 is the random value distributed in [0, 1], c.sub.4 is the acceleration factor related to diversity direction, selected in [0.3, 0.5] randomly, Da(k) is the d-dimensional flight direction of particles with maximum diversity in the population.

[0054] Case 4: When PC(k)<PC(k−1) and PD(k)<PD(k−1), the velocity and position of particle are updated by

v.sub.n,d(k+1)=ωv.sub.n,d(k)+c.sub.1r.sub.1(p.sub.n,d(k)−x.sub.n,d(k))+c.sub.2r.sub.2(g.sub.d(k)−x.sub.i,d(k))+½(c.sub.3r.sub.3C.sub.d(k)+c.sub.4r.sub.dD.sub.d(k))   (20)

x.sub.n,d(k+1)=x.sub.n,d(k)+v.sub.n,d(k+1)  (21)

[0055] Case 5: When PC(k)=PC(k−1) or PD(k)=PD(k−1), the velocity and position of particle are updated by

[00014]$\begin{array}{cc}{v}_{n,d}\left(k+1\right)=\omega v_{n,d}\left(k\right)+c_1{r}_1\left(p_{n,d}\left(k\right)-x_{n,d}\left(k\right)\right)+c_2{r}_2\left(g_d\left(k\right)-x_{i,d}\left(k\right)\right)& \left(22\right)\\ x_{n,d}\left(k+1\right)=\{\begin{array}{cc}{x}_{d,\min }+\left(x_{d,\max }-x_{d,\min }\right)\times U\left(0,1\right),& r_5\le p_b\\ x_{n,d}\left(k\right),& r_5>p_b\end{array}& \left(23\right)\end{array}$

where U(0, 1) is a random value between 0 and 1 which obeys uniform distribution, x.sub.d,min is the minimum boundary value of d-dimensional position, x.sub.d,max is the maximum boundary value of d-dimensional position. When d=1, x.sub.1,min=0.2 mg/L, x.sub.1,max=2 mg/L. When d=2, x.sub.2,min=0.4 mg/L, x.sub.2,max=3 mg/L. r.sub.5 is the random value distributed in [0, 1], p.sub.b is the mutation probability, which is described as

[00015]$\begin{array}{cc}{p}_b=0.5-0.5\times \frac{k}{K}& \left(24\right)\end{array}$

[0056] IV. The population in the kth iteration is combined with the archive A(k−1) to obtain J(k), and then the non-dominated solutions are selected from J(k) to establish A(k);

[0057] V. If k is greater than or equal to K, go to step VI. If k is less than K, go to step I;

[0058] VI. In the archive A(K), a non-dominated solution is randomly selected as the optimal set-point a*(t)=ab(K), ab(K)=[S.sub.NO*(K), S o*(K)], where S.sub.NO*(K) is the optimal set-point of S.sub.NO, S.sub.O*(K) is the optimal set-point of S.sub.O. Then, the optimal set-points are saved;

[0059] (3) Track the optimal set-points S.sub.NO*(K) and S.sub.O*(K):

[0060] PID controller is designed to track the optimal set-points S.sub.NO*(K) and S.sub.O*(K). The expression of PID controller is

[00016]$\begin{array}{cc}\Delta z\left(t\right)=K_p\left[e\left(t\right)+H_l{\int }_0^{t}e\left(t\right)dt+H_d\frac{de\left(t\right)}{dt}\right]& \left(25\right)\\ K_p=\left[\begin{array}{cc}10000& 0\\ 0& 20\end{array}\right]& \left(26\right)\\ H_l=\left[\begin{array}{cc}3000& 0\\ 0& 5\end{array}\right]& \left(27\right)\\ H_d=\left[\begin{array}{cc}100& 0\\ 0& 1\end{array}\right]& \left(28\right)\end{array}$

where Δz(t)=[ΔQ.sub.a(t), ΔK.sub.La5(t)].sup.T is the manipulated variable matrix, ΔQ.sub.a(t) is the change of internal circulation flow in wastewater treatment process, ΔK.sub.La5(t) is the change of oxygen transfer coefficient in the fifth zone of biochemical reactor, K.sub.p is the proportional coefficient matrix, H.sub.l is the integral coefficient matrix, and Ha is the differential coefficient matrix. e(t)=y*(t).sup.T−y(t).sup.T is the control error, y*(t)=[S.sub.NO*(t), S.sub.O*(t)] is the optimal set-point matrix at time t, y(t)=[S.sub.NO(t), S.sub.O(t)] is the actual output matrix.

[0061] (4) The inputs of data-knowledge driven optimal control system of municipal wastewater treatment process are the change of internal circulation flow ΔQ.sub.a(t) and the change of oxygen transfer coefficient in the fifth zone of biochemical reactor ΔK.sub.La5(t). The optimal set-points of S.sub.NO and S.sub.O in municipal wastewater treatment process are tracked and controlled.

[0062] The framework of data-knowledge driven optimal control method is shown in FIG. 1. The tracking result of nitrate nitrogen is shown in FIG. 2. The solid line is the control output and the dotted line is the actual output. X axis shows the time, Y axis shows the concentration of nitrate nitrogen. The tracking error of nitrate nitrogen is shown in FIG. 3. X axis shows the time, Y axis shows the error of nitrate nitrogen. The tracking result of dissolved oxygen is shown in FIG. 4. The solid line is the control output and the dotted line is the actual output. X axis shows the time, Y axis shows the concentration of dissolved oxygen. The tracking error of dissolved oxygen is shown in FIG. 5. X axis shows the time, Y axis shows the error of dissolved oxygen.