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 of a system for wastewater treatment process (WWTP), wherein the system comprises: data acquisition devices, including flow meter, dissolved oxygen detector, nitrate nitrogen detector, ammonia nitrogen detector, suspended solids concentration detector, effluent biochemical oxygen demand detector, effluent chemical oxygen demand detector, effluent ammonia nitrogen detector, effluent total nitrogen detector, effluent suspended solids concentration; control devices, including aeration pump, reflow pump; data transmission devices, including communication interface, remote interface; and a programmable logic controller PLC. the WWTP mainly includes biochemical reaction process and sedimentation process, among them, the biochemical reaction process is divided into five tanks, including an anaerobic tank, an anoxic tank and three aerobic tanks, mainly involving three reaction processes; firstly, wastewater enters an anaerobic tank from a water inlet for reaction of ammonization, secondly, the wastewater enters an anoxic tank from the anaerobic tank for denitrification reaction, finally, the wastewater flows from the anoxic tank to the three aerobic tanks for nitrification reaction, and flows into the secondary sedimentation tank for sedimentation; the anaerobic tank and the anoxic tank are respectively provided with a stirrer to fully mix the activated sludge and the wastewater; the aeration pump arranged at the third aerobic tank can charge oxygen for the three aerobic tanks, and the reflow pump arranged between the third aerobic tank and the anaerobic tank is used for controlling the reflow rate of the internal circulation; the treated wastewater enters a secondary sedimentation tank to realize the separation of activated sludge and effluent; then the effluent is discharged from the upper end of the secondary sedimentation tank, part of the activated sludge is discharged from the bottom end of the secondary sedimentation tank, and the remaining sludge is returned to the anaerobic tank by the control action of the reflow pump arranged at the secondary sedimentation tank; data acquisition sensors are arranged in the anaerobic tank, the anoxic tank, the aerobic tank and the secondary sedimentation tank; the data acquisition device installed in the biochemical reaction process comprises: a flow meter arranged in the anaerobic tank and used for acquiring inflow flow data, a portable dissolved oxygen detector installed in the third aerobic tank for obtaining dissolved oxygen concentration data, a nitrate nitrogen detector arranged in the anoxic tank and used for acquiring nitrate nitrogen concentration data, an ammonia nitrogen detector arranged in the third aerobic tank and used for acquiring ammonia nitrogen data, a suspended matter concentration detector installed in the anoxic tank for acquiring suspended matter concentration data; the data acquisition device installed in the secondary sedimentation tank: an effluent ammonia nitrogen detector for acquiring effluent ammonia nitrogen data, an effluent total nitrogen detector for acquiring effluent total nitrogen data, an effluent suspended matter concentration detector for acquiring effluent suspended matter data, an effluent biochemical oxygen demand detector for obtaining effluent biochemical oxygen demand data, a water outlet chemical oxygen demand detector used for acquiring water outlet chemical oxygen demand data; the data transmission end of each instrument in the wastewater treatment process is connected to the computer by using the communication interface and the remote interface in the data transmission device; the programmable logic controller PLC is used for changing the rotating speed of the driving motor of the aeration pump and the reflow pump to control the output of the aeration pump and the reflow pump; the software module includes data acquisition module, intelligent modeling module, multi-objective optimization module and control module; the data acquisition module is used for acquiring data information such as process variables related to the wastewater treatment process and optimized operation performance indexes; the intelligent modeling module is used for establishing an optimal operation performance index model of the wastewater treatment process, and comprises a pumping energy consumption model, an aeration energy consumption model and an effluent quality model; the multi-objective optimization module is used for simultaneously optimizing a plurality of performance indexes of the wastewater treatment process, including pumping energy, aeration energy and effluent water quality, and obtaining optimal values of related control variables; the control module compares the control action of the aeration pump and the reflow pump with the obtained optimized set value of the control variable, so as to carry out the control action; the dynamic multi-objective particle swarm optimization-based optimal control method for WWTP comprises the following steps: Step 1: acquiring each performance index parameters of WWTP through data acquisition sensors, and transmitting the performance index values and the related process variable values thereof to a data acquisition module, the performance indexes mainly include pumping energy (PE), aeration energy (AE) and effluent quality (EQ); PE mainly refers to the power consumption generated by the internal reflow pump, which is measured by the pumping energy power meter. AE mainly refers to the power consumption generated by aeration pump during aeration, which is measured by aeration energy power meter. EQ mainly refers to the tax or penalty to be paid for discharging pollution to the receiving water body; the pollutants mainly include effluent ammonia nitrogen (S.sub.NH,e), effluent total nitrogen (N.sub.tot,e), effluent suspended solids (SS.sub.e), effluent biochemical oxygen demand (BOD.sub.e), effluent chemical oxygen demand (COD.sub.e) Obtained by the data acquisition devices installed in the secondary sedimentation tank; Step 2: transmit the performance index values and the related process variable values thereof to an intelligent modeling module through a data transmission protocol, and the intelligent modeling module establish a performance index model of the wastewater treatment system according to the related process variables of the performance index obtained by the data acquisition module; the wastewater treatment system performance index model comprises PE model, AE model and EQ model, with the progress of wastewater treatment operation, the collected data is also constantly updated, and performance index models in the intelligent modeling module will be updated and adjusted along with the change of the collected data; the steps for constructing a performance indicator model for sewage treatment process include: {circle around (1)} obtain the operating data of the performance indices and the related process variables, the performance indices mainly contain PE, AE and EQ, the related process variables of PE, AE and EQ contain 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 can be measured by the corresponding sensors; {circle around (2)} establish the nonlinear models of PE, AE, EQ and the related process variables; since the adjusting time of the two controlled variables S.sub.O and S.sub.NO are different, the established models of PE, AE and EQ also have different adjusting periods; because S.sub.O has great influence on AE and EQ, the nonlinear models of AE and EQ will be adjusted every thirty minutes, which is the same as the adjusting period of S.sub.O; in addition, S.sub.NO has great influence on PE, so the model of PE will be adjusted every two hours, which is the same as the adjusting period of S.sub.NO; two sets of objectives are established: $\begin{array}{cc}\{\begin{array}{c}{f}_1\left(t\right)=\underset{r=1}{\overset{10}{.Math.}}{W}_{1r}\left(t\right)\times e^{-{.Math.x\left(t\right)-c_{1r}\left(t\right).Math.}^2/2{b_{1r}\left(t\right)}^2}+W_1\left(t\right)\\ f_2\left(t\right)=\underset{r=1}{\overset{10}{.Math.}}{W}_{2r}\left(t\right)\times e^{-{.Math.l\left(t\right)-c_{2r}\left(t\right).Math.}^2/2{b_{2r}\left(t\right)}^2}+W_2\left(t\right)\end{array}& \left(1\right)\\ \{\begin{array}{c}{f}_1\left(t\right)=\underset{r=1}{\overset{10}{.Math.}}{W}_{1r}\left(t\right)\times e^{-{.Math.x\left(t\right)-c_{1r}\left(t\right).Math.}^2/2{b_{1r}\left(t\right)}^2}+W_1\left(t\right)\\ f_2\left(t\right)=\underset{r=1}{\overset{10}{.Math.}}{W}_{2r}\left(t\right)\times e^{-{.Math.l\left(t\right)-c_{2r}\left(t\right).Math.}^2/2{b_{2r}\left(t\right)}^2}+W_2\left(t\right)\\ f_3\left(t\right)=\underset{r=1}{\overset{10}{.Math.}}{W}_{3r}\left(t\right)\times e^{-{.Math.s\left(t\right)-c_{2r}\left(t\right).Math.}^2/2{b_{3r}\left(t\right)}^2}+W_3\left(t\right)\end{array}& \left(2\right)\end{array}$ Where f.sub.1(t) is AE model at the tth time, f.sub.2(t) is EQ model at the tth time, f.sub.3(t) is PE model at the tth time, x(t)=[Q.sub.in(t), S.sub.O(t)] is the inputs of AE model, l(t)=[Q.sub.in(t), S.sub.O(t), S.sub.NO(t), S.sub.NH(t), SS(t)] is the inputs of EQ model, s(t)=[Q.sub.in(t), S.sub.NO(t), SS(t)] is the inputs of PE model, $e^{-{.Math.x\left(t\right)-c_{1r}\left(t\right).Math.}^2/2{b_{1r}\left(t\right)}^2},e^{-{.Math.l\left(t\right)-c_{2r}\left(t\right).Math.}^2/2{b_{2r}\left(t\right)}^2}\mathrm{and}{e}^{-{.Math.s\left(t\right)-c_{2r}\left(t\right).Math.}^2/2{b_{3r}\left(t\right)}^2}$ are the rth radial basis functions of f.sub.1(t), f.sub.2(t) and f.sub.3(t) at the tth time, r=1, 2, . . . , 10, c.sub.1r(t), c.sub.2r(t) and c.sub.3r(t) are the center vectors of the rth radial basis functions of f.sub.1(t), f.sub.2(t) and f.sub.3(t) at the tth time, and the ranges of c.sub.1r(t), c.sub.2r(t) and c.sub.3r(t) are [−1, 1] respectively; b.sub.1r(t), b.sub.2r(t) and b.sub.3r(t) are the width vectors of the rth radial basis functions of f.sub.1(t), f.sub.2(t) and f.sub.3(t) at the tth time, whose ranges are (0, 2) respectively, W.sub.1r(t), W.sub.2r(t) and W.sub.3r(t) are the weight vectors of the rth radial basis functions of f.sub.1(t) and f.sub.2(t) at the tth time, whose ranges are [−3, 3] respectively; W.sub.1(t), W.sub.2(t) and W.sub.3(t) are the biases of the rth radial basis functions of f.sub.1(t) and f.sub.2(t) at the tth time, whose ranges are [−2, 2] respectively; when every two hours is reached, Eq. (2) is taken as the operating objectives, otherwise, Eq. (1) is regarded as the operating objectives; Step 3: transmitting model results of the PE model, AE model and EQ model built in the intelligent modeling module to the multi-objective optimization module; in the multi-objective optimization module, the constructed PE model, AE model and EQ model are taken as optimization objectives, and then the optimization objectives are optimized by using the designed dynamic multi-objective particle swarm optimization algorithm to obtain the optimized set values of the control variables S.sub.O and S.sub.NO; since the optimal set values of S.sub.O and S.sub.NO change in real time with the progress of wastewater treatment operation process, the optimal set values of S.sub.O and S.sub.NO concentration are recalculated every two hours; the updated optimized set points for dissolved oxygen and nitrate concentration are transmitted to the control module in real time; in WWTP, S.sub.O and S.sub.NO are two important controlled variables to guarantee the treatment results and improve the operating performance; due to the changes of passively accepted influent flow rate, fixed set-points of S.sub.O and S.sub.NO are not beneficial for reducing the PE, AE and EQ. Therefore, a dynamic optimization algorithm is designed for calculating the time-varying set-points of S.sub.O and S.sub.NO; the steps of the dynamic optimization process include: (3)-1 set the maximum iterative numbers of the optimization process T.sub.max; (3)-2 take the established models of performance indices in Eq. (1) and (2) as the optimization objectives; (3)-3 regard all the related inputs of PE, AE and EQ as the position of the particles, a(t)=[Q.sub.in(t), S.sub.O(t), S.sub.NO(t), S.sub.NH(t), SS(t)], 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:
a.sub.k,i(t+1)=a.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 a.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, the range of ω is (0, 1], c.sub.1 and c.sub.2 are the learning parameters, α.sub.1 and α.sub.2 are the uniformly distributed random numbers, 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; (3)-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, $\begin{array}{cc}U\left(t\right)=\sqrt{\frac{1}{\mathrm{NS}\left(t\right)-1}\underset{m=1}{\overset{\mathrm{NS}\left(t\right)}{.Math.}}{\left(\stackrel{⌣}{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 $\begin{array}{cc}A\left(t\right)=\frac{1}{\mathrm{NS}\left(t\right)}\sqrt{\underset{l=1}{\overset{\mathrm{NS}\left(t\right)}{.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 k−1th iteration; (3)-5 judge the changes of the optimization objectives, if the number of the objectives is changed, return to step (3)-6; otherwise, return to (3)-7; (3)-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 $\begin{array}{cc}{N}_{k+1}\left(t\right)=\{\begin{array}{cc}{N}_k\left(t\right)& {\alpha }_k=0\\ N_k\left(t\right)-\left(N_k\left(t\right)-{\mathrm{NS}}_k\left(t\right)\right).Math.{\alpha }_k\left(t\right)& {\alpha }_k<0\\ N_k\left(t\right)+{\mathrm{NS}}_k\left(t\right).Math.{\alpha }_k\left(t\right)& {\alpha }_k>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 $\begin{array}{cc}{\alpha }_k\left(t\right)=\frac{U_k\left(t\right)-U_k\left(t-\epsilon \right)}{\epsilon }& \left(8\right)\end{array}$ where ε is the adjusted frequency of population size, and the 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)& {\beta }_k\left(t\right)=0\\ N_k\left(t\right)+{\mathrm{NS}}_k\left(t\right).Math.{\beta }_k\left(t\right)& {\beta }_k\left(t\right)<0\\ N_k\left(t\right)-\left(N_k\left(t\right)-{\mathrm{NS}}_k\left(t\right)\right).Math.{\beta }_k\left(t\right)& {\beta }_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 $\begin{array}{cc}{\beta }_k\left(t\right)=\frac{A_k\left(t\right)-A_k\left(t-\epsilon \right)}{\epsilon }& \left(10\right)\end{array}$ (3)-7 compare pBest.sub.k(t) with the solutions Φ.sub.k−1(t) in the archive, where Φ.sub.k−1(t)=[φ.sub.k−1,1(t), φ.sub.k−1,2(t), . . . , φ.sub.k−1,i(t)], φ.sub.k−1,i(t) is the tth 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.k−1(t)∪p.sub.k−1(t), if f.sub.h(a.sub.k−1(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.k−1,l(t), then the pBest.sub.k−l(t) will be saved in the archive, otherwise, a.sub.k−1,l(t) will be saved, then gBest.sub.k(t) will be selected from the archive according to the density method; (3)-8 if the current iteration is greater than T.sub.max, then return to step (3)-9, otherwise, return to step (3)-3; (3)-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; Step 4: transmitting the optimized set values of S.sub.O,Tmax* and S.sub.NO,Tmax* obtained in step 3 to a control module, changing the rotating speeds of the driving motors of the aeration pump and the reflow pump by adjusting the variable quantity of the variable frequency output of the programmable logic controller PLC, controlling the air velocity of the aeration pump and the internal circulation flow of the reflow pump, thereby controlling the concentration of S.sub.O in the third aerobic tank and the concentration of S.sub.NO in the anoxic tank. The control action is obtained by a multivariable PID controller, which can be shown as $\begin{array}{cc}\Delta u\left(t\right)=K_p\left[e\left(t\right)+H_{\tau }{\int }_0^{t}e\left(t\right)\mathrm{dt}+H_d\frac{\mathrm{de}\left(t\right)}{\mathrm{dt}}\right]& \left(12\right)\end{array}$ where Δu(t)=[ΔK.sub.La(t), ΔQ.sub.a(t)].sup.T, ΔK.sub.La(t) is the variation of oxygen transfer coefficient in the third aeration tank at time t, ΔQ.sub.a(t) is the variation of 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 errors between the real outputs and the optimal set-points
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) and z.sub.2(t) are the derived optimal set-points of S.sub.O and S.sub.NO at time t, y(t)=[y.sub.1(t), y.sub.2(t)].sup.T, y.sub.1(t) and y.sub.2(t) are the real values of S.sub.O and S.sub.NO at time t; as the optimal set-points of S.sub.O,Tmax* and S.sub.NO,Tmax* are changing timely, it is necessary to compare the control action of airflow rate and internal recycle flow rate with the response of S.sub.O and S.sub.NO concentration; when the real concentration value of S.sub.O is close to the optimized set value of the S.sub.O concentration obtained in the step 3, the control quantity of the current frequency conversion output of the PLC is kept, and the rotating speed of the driving motor of the aeration pump is not changed; when the real value of S.sub.O is less than the optimized set value, the variation of the current frequency conversion output of the PLC is calculated according to the control method in the control module, the rotating speed of the driving motor of the aeration pump is increased, and the aeration amount of the aeration pump is increased, so that the microbial activity is increased and the effluent quality is improved; when the real value of S.sub.O is greater than the optimized set value, the variation of the current frequency conversion output of the PLC is calculated by a control method in the control module, the rotating speed of the driving motor of the aeration pump is reduced, and the aeration amount of the aeration pump is reduced, so that the aeration energy consumption is reduced; when the real concentration value of the nitrate nitrogen is close to the optimized set value of the nitrate nitrogen concentration obtained in the step 3, the control quantity of the current variable frequency output of the PLC is kept, and the rotating speed of the driving motor of the reflow pump is not changed; when the real value of nitrate nitrogen is less than the optimized set value, the variation of the current frequency conversion output of the PLC is calculated according to the variation of the internal circulation flow, the rotating speed of the driving motor of the reflow pump is increased, and the reflow rate of the reflow pump is increased, so that the nitrogen is removed and the effluent quality is improved; when the real value of nitrate nitrogen is greater than the optimized set value, the change amount of the current variable frequency output of the PLC is calculate according to the change of the internal circulation flow, the rotating speed of the driving motor of the reflow pump is reduced, and the reflow rate of the reflow pump is reduced, so that the pumping energy consumption is reduced.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0045] FIG. 1 shows the hardware schematic diagram of the wastewater treatment process, where 1 represents the flow meter, 2 represents a dissolved oxygen detector, 3 represents a nitrate nitrogen detector, 4 represents an ammonia nitrogen detector, 5 represents a suspended solids concentration detector, 6 refers to effluent biochemical oxygen demand detector, 7 refers to the chemical oxygen demand detector of effluent, 8 represents the effluent ammonia nitrogen detector, 9 represents the total nitrogen detector for effluent, 10 represents the concentration of suspended solids in the effluent, 11 represents the pumping energy consumption power meter, 12 represents the aeration energy consumption power meter.

[0046] FIG. 2 shows the results of S.sub.O in DMOPSO-based method.

[0047] FIG. 3 shows the results of S.sub.NO in DMOPSO-based method.

DETAILED DESCRIPTION

[0048] In this present invention, a dynamic multi-objective particle swarm optimization-based optimal control method and system for wastewater treatment process is designed, where hardware facilities and software modules are included.

[0049] The hardware facilities include data acquisition devices (flow meter, dissolved oxygen detector, nitrate nitrogen detector, ammonia nitrogen detector, suspended solids concentration detector, effluent biochemical oxygen demand detector, effluent chemical oxygen demand detector, effluent ammonia nitrogen detector, effluent total nitrogen detector, effluent suspended solids concentration), control devices (aeration pump, reflow pump), data transmission devices (communication interface, remote interface) and programmable logic controller PLC.

[0050] As the controlled object, WWTP mainly includes biochemical reaction process and sedimentation process, among them, the biochemical reaction process is divided into five tanks, including an anaerobic tank, an anoxic tank and three aerobic tanks, mainly involving three reaction processes. Firstly, wastewater enters an anaerobic tank from a water inlet for reaction of ammonization. Secondly, the wastewater enters an anoxic tank from the anaerobic tank for denitrification reaction. Finally, the wastewater flows from the anoxic tank to the three aerobic tanks for nitrification reaction, and flows into the secondary sedimentation tank for sedimentation. The anaerobic tank and the anoxic tank are respectively provided with a stirrer to fully mix the activated sludge and the wastewater. The aeration pump arranged at the third aerobic tank can charge oxygen for the three aerobic tanks, and the reflow pump arranged between the third aerobic tank and the anaerobic tank is used for controlling the reflow rate of the internal circulation. The treated wastewater enters a secondary sedimentation tank to realize the separation of activated sludge and effluent. Then the effluent is discharged from the upper end of the secondary sedimentation tank, part of the activated sludge is discharged from the bottom end of the secondary sedimentation tank, and the remaining sludge is returned to the anaerobic tank by the control action of the reflow pump arranged at the secondary sedimentation tank. Data acquisition sensors are arranged in the anaerobic tank, the anoxic tank, the aerobic tank and the secondary sedimentation tank.

[0051] The data acquisition device installed in the biochemical reaction process comprises:

[0052] A flow meter arranged in the anaerobic tank and used for acquiring inflow flow data, a portable dissolved oxygen detector installed in the third aerobic tank for obtaining dissolved oxygen concentration data, a nitrate nitrogen detector arranged in the anoxic tank and used for acquiring nitrate nitrogen concentration data, an ammonia nitrogen detector arranged in the third aerobic tank and used for acquiring ammonia nitrogen data, a suspended matter concentration detector installed in the anoxic tank for acquiring suspended matter concentration data.

[0053] The data acquisition device installed in the secondary sedimentation tank:

[0054] An effluent ammonia nitrogen detector for acquiring effluent ammonia nitrogen data, an effluent total nitrogen detector for acquiring effluent total nitrogen data, an effluent suspended matter concentration detector for acquiring effluent suspended matter data, an effluent biochemical oxygen demand detector for obtaining effluent biochemical oxygen demand data, a water outlet chemical oxygen demand detector used for acquiring water outlet chemical oxygen demand data.

[0055] The data transmission end of each instrument in the wastewater treatment process is connected to the computer by using the communication interface and the remote interface in the data transmission device. The programmable logic controller PLC is used for changing the rotating speed of the driving motor of the aeration pump and the reflow pump to control the output of the aeration pump and the reflow pump.

[0056] The software module includes data acquisition module, intelligent modeling module, multi-objective optimization module and control module.

[0057] The data acquisition module is used for acquiring data information such as process variables related to the wastewater treatment process and optimized operation performance indexes.

[0058] The intelligent modeling module is used for establishing an optimal operation performance index model of the wastewater treatment process, and comprises a pumping energy consumption model, an aeration energy consumption model and an effluent quality model.

[0059] The multi-objective optimization module is used for simultaneously optimizing a plurality of performance indexes of the wastewater treatment process, including pumping energy, aeration energy and effluent water quality, and obtaining optimal values of related control variables.

[0060] The control module compares the control action of the aeration pump and the reflow pump with the obtained optimized set value of the control variable, so as to carry out the control action.

[0061] The invention relates to a dynamic multi-objective particle swarm optimization-based optimal control method for WWTP, which comprises the following steps:

[0062] Step 1: acquiring each performance index parameters of WWTP through data acquisition sensors, and transmitting the performance index values and the related process variable values thereof to a data acquisition module. The performance indexes mainly include pumping energy (PE), aeration energy (AE) and effluent quality (EQ).

[0063] PE mainly refers to the power consumption generated by the internal reflow pump, which is measured by the pumping energy power meter. AE mainly refers to the power consumption generated by aeration pump during aeration, which is measured by aeration energy power meter. EQ mainly refers to the tax or penalty to be paid for discharging pollution to the receiving water body. The pollutants mainly include effluent ammonia nitrogen (S.sub.NH,e), effluent total nitrogen (N.sub.tot,e), effluent suspended solids (SS.sub.e), effluent biochemical oxygen demand (BOD.sub.e), effluent chemical oxygen demand (COD.sub.e) obtained by the data acquisition devices installed in the secondary sedimentation tank.

[0064] Step 2: transmit the performance index values and the related process variable values thereof to an intelligent modeling module through a data transmission protocol, and the intelligent modeling module establish a performance index model of the wastewater treatment system according to the related process variables of the performance index obtained by the data acquisition module.

[0065] The wastewater treatment system performance index model comprises PE model, AE model and EQ model. With the progress of wastewater treatment operation, the collected data is also constantly updated, and performance index models in the intelligent modeling module will be updated and adjusted along with the change of the collected data. The steps for constructing a performance indicator model for sewage treatment process include:

[0066] {circle around (1)} obtain the operating data of the performance indices and the related process variables. The performance indices mainly contain PE, AE and EQ. The related process variables of PE, AE and EQ contain 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 can be measured by the corresponding sensors.

[0067] {circle around (2)} establish the nonlinear models of PE, AE, EQ and the related process variables. Since the adjusting time of the two controlled variables S.sub.O and S.sub.NO are different, the established models of PE, AE and EQ also have different adjusting periods. Because S.sub.O has great influence on AE and EQ, the nonlinear models of AE and EQ will be adjusted every thirty minutes, which is the same as the adjusting period of S.sub.O. In addition, S.sub.NO has great influence on PE, so the model of PE will be adjusted every two hours, which is the same as the adjusting period of S.sub.NO. Two sets of objectives are established

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

[0068] where f.sub.1(t) is AE model at the tth time, f.sub.2(t) is EQ model at the tth time, f.sub.3(t) is PE model at the tth time, x(t)=[Q.sub.in(t), S.sub.O(t)] is the inputs of AE model, l(t)=[Q.sub.in(t), S.sub.O(t), S.sub.NO(t), S.sub.NH(t), SS(t)] is the inputs of EQ model, s(t)=[Q.sub.in(t), S.sub.NO(t), SS(t)] is the inputs of PE model,

[00011]$e^{-{.Math.x\left(t\right)-c_{1r}\left(t\right).Math.}^2/2{b_{1r}\left(t\right)}^2},e^{-{.Math.l\left(t\right)-c_{2r}\left(t\right).Math.}^2/2{b_{2r}\left(t\right)}^2}\mathrm{and}{e}^{-{.Math.s\left(t\right)-c_{3r}\left(t\right).Math.}^2/2{b_{3r}\left(t\right)}^2}$

are the rth radial basis functions of f.sub.1(t), f.sub.2(t) and f.sub.3(t) at the tth time, r=1, 2, . . . , 10, c.sub.1r(t), c.sub.2r(t) and c.sub.3r(t) are the center vectors of the rth radial basis functions of f.sub.1(t), f.sub.2(t) and f.sub.3(t) at the tth time, and the ranges of c.sub.1r(t), c.sub.2r(t) and c.sub.3r(t) are [−1, 1] respectively; b.sub.1r(t), b.sub.2r(t) and b.sub.3r(t) are the width vectors of the rth radial basis functions of f.sub.1(t), f.sub.2(t) and f.sub.3(t) at the tth time, whose ranges are (0, 2] respectively, W.sub.1r(t), W.sub.2r(t) and W.sub.3r(t) are the weight vectors of the rth radial basis functions of f.sub.1(t) and f.sub.2(t) at the tth time, whose ranges are [−3, 3] respectively; W.sub.1(t), W.sub.2(t) and W.sub.3(t) are the biases of the rth radial basis functions of f.sub.1(t) and f.sub.2(t) at the tth time, whose ranges are [−2, 2] respectively. When every two hours is reached, Eq. (2) is taken as the operating objectives, otherwise, Eq. (1) is regarded as the operating objectives.

[0069] Step 3: transmitting model results of the PE model, AE model and EQ model built in the intelligent modeling module to the multi-objective optimization module. In the multi-objective optimization module, the constructed PE model, AE model and EQ model are taken as optimization objectives, and then the optimization objectives are optimized by using the designed dynamic multi-objective particle swarm optimization algorithm to obtain the optimized set values of the control variables S.sub.O and S.sub.NO.

[0070] Since the optimal set values of S.sub.O and S.sub.NO change in real time with the progress of wastewater treatment operation process, the optimal set values of S.sub.O and S.sub.NO concentration are recalculated every two hours. The updated optimized set points for dissolved oxygen and nitrate concentration are transmitted to the control module in real time.

[0071] In WWTP, S.sub.O and S.sub.NO are two important controlled variables to guarantee the treatment results and improve the operating performance. Due to the changes of passively accepted influent flow rate, fixed set-points of S.sub.O and S.sub.NO are not beneficial for reducing the PE, AE and EQ. Therefore, a dynamic optimization algorithm is designed for calculating the time-varying set-points of S.sub.O and S.sub.NO. The steps of the dynamic optimization process include:

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

[0073] (3)-2 take the established models of performance indices in Eq. (1) and (2) as the optimization objectives;

[0074] (3)-3 regard all the related inputs of PE, AE and EQ as the position of the particles, a(t)=[Q.sub.in(t), S.sub.O(t), S.sub.NO(t), S.sub.NH(t), SS(t)], 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:

a.sub.k,i(t+1)=a.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 a.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;

[0075] (3)-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}\underset{m=1}{\overset{\mathrm{NS}\left(t\right)}{.Math.}}{\left(\stackrel{⌣}{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, {circle around (D)}(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)}\sqrt{\underset{l=1}{\overset{\mathrm{NS}\left(t\right)}{.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 k−1th iteration;

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

[0077] (3)-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)& {\alpha }_k=0\\ N_k\left(t\right)-\left(N_k\left(t\right)-{\mathrm{NS}}_k\left(t\right)\right).Math.{\alpha }_k\left(t\right)& {\alpha }_k<0\\ N_k\left(t\right)+{\mathrm{NS}}_k\left(t\right).Math.{\alpha }_k\left(t\right)& {\alpha }_k>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, a.sub.k(t) is the gradient of diversity in the kth iteration at the tth time, which is calculated as

[00015]$\begin{array}{cc}{\alpha }_k\left(t\right)=\frac{U_k\left(t\right)-U_k\left(t-\epsilon \right)}{\epsilon }& \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)& {\beta }_k\left(t\right)=0\\ N_k\left(t\right)+{\mathrm{NS}}_k\left(t\right).Math.{\beta }_k\left(t\right)& {\beta }_k\left(t\right)<0\\ N_k\left(t\right)-\left(N_k\left(t\right)-{\mathrm{NS}}_k\left(t\right)\right).Math.{\beta }_k\left(t\right)& {\beta }_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}{\beta }_k\left(t\right)=\frac{A_k\left(t\right)-A_k\left(t-\epsilon \right)}{\epsilon }& \left(10\right)\end{array}$ [0078] (3)-7 compare pBest.sub.k(t) with the solutions Φ.sub.k−1(t) in the archive, where Φ.sub.k−1(t)=[φ.sub.k−1,1(t), φ.sub.k−1,2(t), . . . , φ.sub.k−i,l(t)], φ.sub.k−1,i(t) is the lth 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.k−1(t)∪p.sub.k−1(t), if f.sub.h(a.sub.k−1(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.k−1,l(t), then the pBest.sub.k−1(t) will be saved in the archive, otherwise, a.sub.k−1,l(t) will be saved, then gBest.sub.k(t) will be selected from the archive according to the density method;

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

[0080] (3)-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.

[0081] Step 4: transmitting the optimized set values of S.sub.O,Tmax* and S.sub.NO,Tmax* obtained in step 3 to a control module, changing the rotating speeds of the driving motors of the aeration pump and the reflow pump by adjusting the variable quantity of the variable frequency output of the programmable logic controller PLC, controlling the air velocity of the aeration pump and the internal circulation flow of the reflow pump, thereby controlling the concentration of S.sub.O in the third aerobic tank and the concentration of S.sub.NO in the anoxic tank. The control action is obtained by a multivariable PID controller, which can be shown as

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

where Δu(t)=[ΔKLa(t), ΔQ.sub.a(t)].sup.T, ΔK.sub.La(t) is the variation of oxygen transfer coefficient in the third aeration tank at time t, ΔQ.sub.a(t) is the variation of 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) and z.sub.2(t) are the derived optimal set-points of S.sub.O and S.sub.NO at time t, y(t)=[y.sub.1(t), y.sub.2(t)].sup.T, y.sub.1(t) and y.sub.2(t) are the real values of S.sub.O and S.sub.NO at time t.

[0082] As the optimal set-points of S.sub.O,Tmax* and S.sub.NO,Tmax* are changing timely, it is necessary to compare the control action of airflow rate and internal recycle flow rate with the response of S.sub.O and S.sub.NO concentration.

[0083] When the real concentration value of S.sub.O is close to the optimized set value of the S.sub.O concentration obtained in the step 3, the control quantity of the current frequency conversion output of the PLC is kept, and the rotating speed of the driving motor of the aeration pump is not changed. When the real value of S.sub.O is less than the optimized set value, the variation of the current frequency conversion output of the PLC is calculated according to the control method in the control module, the rotating speed of the driving motor of the aeration pump is increased, and the aeration amount of the aeration pump is increased, so that the microbial activity is increased and the effluent quality is improved. When the real value of S.sub.O is greater than the optimized set value, the variation of the current frequency conversion output of the PLC is calculated by a control method in the control module, the rotating speed of the driving motor of the aeration pump is reduced, and the aeration amount of the aeration pump is reduced, so that the aeration energy consumption is reduced.

[0084] When the real concentration value of the nitrate nitrogen is close to the optimized set value of the nitrate nitrogen concentration obtained in the step 3, the control quantity of the current variable frequency output of the PLC is kept, and the rotating speed of the driving motor of the reflow pump is not changed. When the real value of nitrate nitrogen is less than the optimized set value, the variation of the current frequency conversion output of the PLC is calculated according to the variation of the internal circulation flow, the rotating speed of the driving motor of the reflow pump is increased, and the reflow rate of the reflow pump is increased, so that the nitrogen is removed and the effluent quality is improved. When the real value of nitrate nitrogen is greater than the optimized set value, the change amount of the current variable frequency output of the PLC is calculate according to the change of the internal circulation flow, the rotating speed of the driving motor of the reflow pump is reduced, and the reflow rate of the reflow pump is reduced, so that the pumping energy consumption is reduced.

[0085] The outputs of the designed DMOPSO-based optimal control method is the concentration of S.sub.O and S.sub.NO. FIG. 2(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. 2(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. 3(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. 3(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.