Urban sewage treatment process time delay robust optimization control method
Patent Information
- Application Number
- CN202310614392.X
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- Filing Date
- 2023-05-29
- Publication Date
- 2026-10-09
- Estimated Expiration
- 2043-05-29
AI Technical Summary
然而,由于反应物浓度变化在污水处理过程中存在时间滞后性,使得对其的测量存在偏差,进而造成了优化目标模型的不确定性,难以实现精确的在线优化控制,常常发生出水水质超标和运行能耗过高的现象,给城市环境和污水处理厂都带来了较重的负担
Smart Images

Figure CN116776722B_ABST
Abstract
Description
Technical Field
[0001] This invention addresses the time-delay characteristics of urban wastewater treatment processes by proposing a time-delay robust optimization control method. This method overcomes the influence of time-delay deviations in biochemical reactant concentration measurements, enhances the robustness of optimized setpoints for dissolved oxygen and nitrate nitrogen, and achieves more precise tracking control. It reduces operating energy consumption while ensuring effluent quality. Controlling effluent quality and operating energy consumption are core objectives of urban wastewater treatment plant operation optimization and are of great significance for the efficient and stable operation of urban wastewater treatment processes. This method belongs to both the field of intelligent optimization control and the field of water treatment. Background Technology
[0002] With the continuous increase in urban population, the demand for freshwater resources is constantly rising, and the problem of urban sewage treatment is becoming increasingly prominent. To address water scarcity and water pollution, the State Council has issued documents such as the "Action Plan for Water Pollution Prevention and Control," encouraging and supporting scientific and technological research on water pollution prevention and control, as well as the promotion and application of advanced and applicable technologies. Urban sewage treatment can purify water quality, improve the ecological environment, and promote water resource recycling, which is of great significance to urban development. Improving effluent quality and reducing operating energy consumption are the core objectives of optimizing the operation of urban sewage treatment plants, and have a significant impact on protecting the urban environment and enhancing the social benefits of sewage treatment plants. Therefore, research on the optimal control of effluent quality and operating energy consumption in urban sewage treatment processes has broad application prospects.
[0003] With the continuous improvement of urban wastewater treatment optimization and control technologies, organic pollutants in wastewater have been effectively controlled. However, due to the time lag in reactant concentration changes during wastewater treatment, measurement deviations occur, leading to uncertainties in the optimization target model. This makes precise online optimization and control difficult, often resulting in effluent quality exceeding standards and excessive energy consumption, placing a heavy burden on the urban environment and wastewater treatment plants. Therefore, enhancing the robustness of dissolved oxygen and nitrate nitrogen optimization setpoints under the influence of the time lag characteristics of wastewater treatment, and stably ensuring effluent quality while reducing operational energy consumption, is an important research topic.
[0004] This invention designs a time-delay robust optimization control method for urban wastewater treatment processes. It mainly analyzes the characteristics of the urban wastewater treatment process, establishes a time-delay multi-objective optimization model of effluent quality and operating energy consumption, and designs a robust evolutionary algorithm to solve the model. This yields robust desired setpoints for dissolved oxygen and nitrate nitrogen, and controls the process variables to track these setpoints. This overcomes the influence of time-delay disturbances during wastewater treatment, reduces operating energy consumption while ensuring effluent quality, and achieves efficient and stable operation of the wastewater treatment process. Summary of the Invention
[0005] This invention designs a time-delay robust optimization control method for urban wastewater treatment processes. This method combines data-driven modeling and time-delay compensation methods to establish a time-delay multi-objective optimization model for effluent quality and operating energy consumption in urban wastewater treatment processes. By designing a time-delay robust evolutionary optimization algorithm, robust setpoints for dissolved oxygen and nitrate nitrogen are solved. A proportional-integral-derivative controller is used to track the expected optimal setpoints of process variables in real time. This overcomes the uncertainties caused by time-delay disturbances at the biochemical reaction end and the effluent end of the wastewater treatment process, effectively improving the effluent quality and reducing operating energy consumption.
[0006] The present invention adopts the following technical solution and implementation steps:
[0007] 1. A time-delay robust optimization control method for urban wastewater treatment processes, specifically including the following steps:
[0008] (1) Data collection of key process variables: Seven variables were selected as key process variables: dissolved oxygen concentration in the aerobic reactor, nitrate nitrogen concentration in the anoxic reactor, ammonia nitrogen concentration in the effluent, total nitrogen concentration in the effluent, total suspended solids concentration in the effluent, bio-oxygen demand concentration in the effluent, and chemical oxygen demand concentration in the effluent; the key process variable data were collected by sensors at time t as follows:
[0009] S(t)=[S O (t),S NO (t),S eNH (t),S eTN (t),S eTSS (t),S eBOD (t),S eCOD (t)] (1)
[0010] Where S(t) is the key variable vector at time t, S O (t) represents the dissolved oxygen concentration in the aerobic reactor at time t, in mg / L, S. NO (t) represents the nitrate nitrogen concentration in the anoxic reactor at time t, in mg / L, S eNH (t) represents the ammonia nitrogen concentration in the effluent at time t, in mg / L, S eTN (t) represents the total nitrogen concentration in the effluent at time t, in mg / L, S eTSS (t) represents the total suspended solids concentration in the effluent at time t, in mg / L, S eBOD (t) represents the biological oxygen demand at time t, in mg / L, S eCOD (t) represents the chemical oxygen demand at time t, in mg / L;
[0011] (2) Establish a time-delay robust multi-objective optimization model for urban wastewater treatment process: Use kernel functions to establish a time-delay robust multi-objective optimization model for urban wastewater treatment process, including effluent quality optimization model, operation energy consumption optimization model and time-delay robust compensation model, to express the relationship between key process variables and effluent quality and operation energy consumption;
[0012] ① Establish an effluent water quality optimization model:
[0013]
[0014] Among them, J EQ (x(t)) is the effluent water quality optimization model at time t, x(t) = [S O (t),S NO [t] represents J at time t. EQ The input vector of (x(t)), δ EQ,k (t) represents J at time t. EQ The k-th weight of (x(t)) and its value range is [-10, 10], c EQ,k (t)=[c EQ,k,1 (t),c EQ,k,2 [t] represents J at time t. EQ The k-th center vector of (x(t)), c EQ,k,1 (t) is c EQ,k The first component of (t) takes values in the range [-5, 5], c EQ,k,2 (t) is c EQ,k The second component of (t) takes values in the range [-5, 5], b EQ,k (t) represents J at time t. EQ The k-th width of (x(t)) and its value range is [0,3];
[0015] ② Establish an operational energy consumption optimization model:
[0016]
[0017] Among them, J EC (x(t)) is the energy consumption optimization model at time t, x(t) = [S O (t),S NO [t] represents J at time t. EC The input vector of (x(t)), δ EC,k (t) represents J at time t. EC The k-th weight of (x(t)) and its value range is [-10, 10], c EC,k (t)=[c EC,k,1 (t),c EC,k,2 [t] represents J at time t. EC The k-th center vector of (x(t)), cEC,k,1 (t) is c EC,k The first component of (t) takes values in the range [-5, 5], c EC,k,2 (t) is c EC,k The second component of (t) takes values in the range [-5, 5], b EC,k (t) represents J at time t. EC The k-th width of (x(t)) and its value range is [0,3];
[0018] ③ Establish a time-delay robust compensation model:
[0019]
[0020] Wherein, ε(x) τ (t) represents the time-delay robust compensation model at time t, x τ (t)=[S O (t-τ DO ),S NO (t-τ NO )] represents ε(x) at time t τ The input vector of τ(t)) DO =0.5(τ) DO,low +τ DO,up τ represents the time lag of dissolved oxygen in the aerobic reactor relative to the effluent. DO,low =V S / Q in For τ DO The lower bound, V S Q is the volume of the secondary sedimentation tank. in τ is the volume of water entering the vessel. DO,up =(V S +V O ) / Q in For τ DO The upper bound of V O S is the volume of the aerobic reactor. O (t-τ DO ) is t-τ DO The dissolved oxygen concentration in the aerobic reactor at any given time is within the range of [0,3], and the unit is mg / L; S NO (t-τ NO ) is t-τ NO The nitrate nitrogen concentration in the aerobic reactor at a given time is within the range of [0,2], and the unit is mg / L. NO =0.5(τ) DO,low +τ DO,up τ represents the time lag of nitrate nitrogen relative to the effluent in the anoxic reactor. NO,low =(V S +3V O ) / Q in For τNO The lower bound of τ NO,up =(V S +V A +3V O ) / Q in For τ NO The upper bound of V A The volume of the anoxic reactor; δ ε,k (t) represents ε(x) at time t. τ The k-th weight of (t) is in the range [-10, 10], c ε,k (t)=[c ε,k,1 (t),c ε,k,2 [(t)] represents ε(x) at time t. τ The k-th center vector of (t) is c ε,k,1 (t) is c ε,k The first component of (t) takes values in the range [-5, 5], c ε,k,2 (t) is c ε,k The second component of (t) takes values in the range [-5, 5], b ε,k (t) represents ε(x) at time t. τ The k-th width of (t) is in the range of [0,3].
[0021] ④ For J EQ (x(t)), J EC (x(t)) and ε(x) τ The parameters in (t) are updated in real time:
[0022]
[0023]
[0024]
[0025]
[0026]
[0027]
[0028]
[0029]
[0030]
[0031] Where, δ EQ (t)=[δ EQ,1 (t),δ EQ,2 (t),…,δEQ,5 [t] represents J at time t. EQ The weight vector of (x(t)), δ EC (t)=[δ EC,1 (t),δ EC,2 (t),…,δ EC,5 [t] represents J at time t. EC The weight vector of (x(t)), δ ε (t)=[δ ε,1 (t),δ ε,2 (t),…,δ ε,5 [(t)] represents ε(x) at time t. τ The weight vector of (t)), c EQ (t)=[c EQ,1 (t),c EQ,2 (t),…,c EQ,5 [t] represents J at time t. EQ The center vector of (x(t)), c EC (t)=[c EC,1 (t),c EC,2 (t),…,c EC,5 [t] represents J at time t. EC The center vector of (x(t)), c ε (t)=[c ε,1 (t),c ε,2 (t),…,c ε,5 [(t)] represents ε(x) at time t. τ The center vector of (t)), b EQ (t)=[b EQ,1 (t),b EQ,2 (t),…,b EQ,5 [t] represents J at time t. EQ The width vector of (x(t)), b EC (t)=[b EC,1 (t),b EC,2 (t),…,b EC,5 [t] is J EC The width vector of (x(t)) at time t, b ε (t)=[b ε,1 (t),b ε,2 (t),…,b ε,5 [(t)] represents ε(x) at time t. τ The width vector of (t)), δ EQ (t+1)=[δ EQ,1 (t+1),δ EQ,2 (t+1),…,δ EQ,5 [t+1] represents J at time t+1. EQThe weight vector of (x(t+1)), δ EC (t+1)=[δ EC,1 (t+1),δ EC,2 (t+1),…,δ EC,5 [t+1] represents J at time t+1. EC The weight vector of (x(t+1)), δ ε (t+1)=[δ ε,1 (t+1),δ ε,2 (t+1),…,δ ε,5 [(t+1)] represents ε(x) at time t+1. τ The weight vector of (t+1)), c EQ (t+1)=[c EQ,1 (t+1),c EQ,2 (t+1),…,c EQ,5 [t+1] represents J at time t+1. EQ The center vector of (x(t+1)), c EC (t+1)=[c EC,1 (t+1),c EC,2 (t+1),…,c EC,5 [t+1] represents J at time t+1. EC The center vector of (x(t+1)), c ε (t+1)=[c ε,1 (t+1),c ε,2 (t+1),…,c ε,5 [(t+1)] represents ε(x) at time t+1. τ The center vector of (t+1)), b EQ (t+1)=[b EQ,1 (t+1),b EQ,2 (t+1),…,b EQ,5 [t+1] represents J at time t+1. EQ The width vector of (x(t+1)), b EC (t+1)=[b EC,1 (t+1),b EC,2 (t+1),…,b EC,5 [t+1] represents J at time t+1. EC The width vector of (x(t+1)), b ε (t+1)=[b ε,1 (t+1),b ε,2 (t+1),…,b ε,5 [(t+1)] represents ε(x) at time t+1. τ (t+1)) width vector; e EQ (t)=J EQ(x(t))-EQ(t) represents the training error of the effluent water quality optimization model at time t, EQ(t) = S eNH (t)+S eTN (t)+S eTSS (t)+S eBOD (t)+S eCOD (t) represents the effluent water quality sample value at time t, e EC (t)=J EC (x(t))-EC(t) represents the training error of the energy consumption optimization model at time t, and EC(t) represents the energy consumption collected at time t. ε (t)=e EQ (t)-ε(x τ (t) represents the time-delay robust compensation model training error at time t;
[0032] ⑤ Establish a time-delay robust multi-objective optimization function for the urban wastewater treatment process:
[0033]
[0034] J EQ (x(t)) represents the objective function value of the decision variable x(t) at time t, J EC (x(t)) represents the objective function value of the decision variable x(t) at time t, ε(x τ (t) represents the decision variable x. τ (t) The time delay robust compensation value at time t;
[0035] (3) Solving for optimal setpoints of process variables: Designing a time-delay robust evolutionary optimization algorithm to solve for optimal setpoints of dissolved oxygen and nitrate nitrogen;
[0036] ① Set the total number of iterations for the time-delay robust evolutionary optimization algorithm to κ. max =500, the particle swarm size is N=100, and the particle positions and velocities are:
[0037] z t,n (κ)=[z t,n,1 (κ),z t,n,2 (κ)](15)
[0038] v t,n (κ)=[v t,n,1 (κ),v t,n,2 (κ)](16) where z t,n (κ) is the position vector of the nth particle at time t, which evolves to the κth generation. t,n,1 (κ) is z t,n The first component of (κ) represents the dissolved oxygen concentration, and its value ranges from [0,3]. t,n,2(κ) is z t,n The second component of (κ) represents the nitrate nitrogen concentration, with a value range of [0,2]. t,n (κ) is the velocity vector of the nth particle evolving to the κth generation at time t, v t,n,1 (κ) is v t,n The first component of (κ) and its value range is [0, 0.1], v t,n,2 (κ) is v t,n The second component of (κ) and its value range is [0, 0.1];
[0039] ② Calculate the individual optimal position of each particle and the global optimal position of the population:
[0040]
[0041]
[0042] Where, p t,n (κ+1) represents the optimal position of the nth particle at time t in the (κ+1)th generation, p t,n (κ) represents the optimal position of the nth particle at time t in the historical individual evolution to the κth generation, z t,n (κ+1) is the position vector of the nth particle evolving to the (κ+1)th generation at time t, and g t (κ+1) represents the globally optimal position of the population at time t after evolving to the (κ+1)th generation, and g t (κ) represents the globally optimal position of the population at time t in the κth generation, and > represents a robust dominance relationship defined as:
[0043]
[0044]
[0045] Among them, J EQ (z t,n (κ+1) is the decision variable z. t,n The objective function value for optimizing effluent water quality is (κ+1); J EC (z t,n (κ+1) is the decision variable z. t,n The objective function value for optimizing the operating energy consumption of (κ+1); ε(z) t,n (κ+1) is the decision variable z. t,n (κ+1) time-delay robust compensation value; J EQ (p t,n (κ) is the decision variable p t,n The objective function value for optimizing effluent water quality (κ); J EC (p t,n (κ) is the decision variable p t,nThe objective function value for optimizing the operating energy consumption of (κ); ε(p t,n (κ)) is the decision variable (p) t,n (κ)) time delay robust compensation value; J EQ (g t,n (κ) is the decision variable g t,n The objective function value for optimizing effluent water quality (κ); J EC (g t,n (κ) is the decision variable g t,n The objective function value for optimizing the operating energy consumption of (κ); ε(g) t,n (κ)) is the decision variable (g) t,n (κ)) time delay robust compensation value; The expected value of the time-delay robust compensation model in the robust region:
[0046]
[0047] Ω=(S NO,up (t)-S NO,low (t))(S O,up (t)-S O,low (t)) (22)
[0048]
[0049]
[0050]
[0051]
[0052] Wherein, ε(z) t,n (κ) is the decision variable z t,n The time-delay robust compensation value of (κ), S NO,up (t) represents S at time t. NO The upper bound of (t), S NO,low (t) represents S at time t. NO The lower bound of S(t) O,up (t) represents S at time t. O The upper bound of (t), S O,low (t) represents S at time t. O (t) is the lower bound of , where Ω is the area of the robust region and λ is the sampling time;
[0053] ③ Update the particle's velocity and position:
[0054] z t,n (κ+1)=z t,n (κ)+v t,n (κ+1) (27)
[0055] v t,n (κ+1)=0.7v t,n (κ)+0.5μ1(p t,n (κ)-v t,n (κ))+0.5μ2(g t (κ)-v t,n (κ)) (28)
[0056] Among them, v t,n (κ+1) is the velocity vector of the nth particle evolving to the κ+1th generation at time t, μ1 is the random number of individual evolution with a value range of [0,1], and μ2 is the random number of population evolution with a value range of [0,1].
[0057] ④ Check if the evolutionary process has reached the stopping condition: If the evolutionary generation κ = 500, terminate the evolutionary process and jump to step ⑤; otherwise, increment the evolutionary generation κ by 1 and return to step ②.
[0058] ⑤ From the κth max The solution set Z of the generation t (κ max )=[z t,1 (κ max ),z t,2 (κ max ),…,z t,N (κ max A solution x is randomly selected from [ ] ] * (t)=[S * O (t),S * NO [(t)] serves as the optimized setpoint for dissolved oxygen and nitrate nitrogen, where S * O (t) represents the optimal setpoint for dissolved oxygen in the aerobic reactor at time t, S * NO (t) represents the optimal setpoint for nitrate nitrogen in the anoxic reactor at time t;
[0059] (4) Tracking control optimization setpoints:
[0060] ① Solve for the optimal setpoint x using a multivariable proportional-integral-derivative control algorithm. * (t) Perform tracking control:
[0061]
[0062] Where, Δu(t)=[ΔK L a(t), ΔQ a (t)] T Let ΔK be the matrix of operations at time t. La(t) represents the change in the dissolved oxygen transfer coefficient at time t, ΔQ a (t) represents the change in circulating flow rate at time t, K p H is a proportionality constant with a value range of [0,5]. l H is the integral coefficient and its value ranges from [0,5]. d The differential coefficients are in the range [0,5], e(t) = x * (t)-x(t) is the control error vector at time t;
[0063] ② Adjust the dissolved oxygen transfer coefficient and internal reflux flow rate using Δu(t):
[0064] K L a(t+1)=K L a(t)+ΔK L a(t) (30)
[0065] Q a (t+1)=Q a (t)+ΔQ a (t) (31)
[0066] Among them, K L a(t+1) is the dissolved oxygen transfer coefficient at time t+1, K L a(t) is the dissolved oxygen transfer coefficient at time t, Q a (t+1) represents the internal return flow rate at time t+1, Q a (t) represents the internal return flow rate at time t;
[0067] (5) K is transmitted via the host computer L a(t) and Q a The value of (t) is transmitted to the programmable logic controller (PLC), and the digital signal is transmitted to the frequency converter via a digital-to-analog converter (DAC). The frequency converter adjusts the frequency of the oxygen supply pump and the internal return pump, thereby adjusting the dissolved oxygen concentration to S. * O (t), the nitrate nitrogen concentration was adjusted to S * NO (t). Attached Figure Description
[0068] Figure 1 The present invention relates to the effluent quality of the urban wastewater treatment process;
[0069] Figure 2 This invention relates to the energy consumption of the urban wastewater treatment process;
[0070] Figure 3 These are the dissolved oxygen concentration control effect diagram and error diagram of the present invention;
[0071] Figure 4These are the effect diagram and error diagram of the nitrate nitrogen concentration control of the present invention. Detailed Implementation
[0072] 1. A time-delay robust optimization control method for urban wastewater treatment processes, specifically including the following steps:
[0073] (1) Data collection of key process variables: Seven variables were selected as key process variables: dissolved oxygen concentration in the aerobic reactor, nitrate nitrogen concentration in the anoxic reactor, ammonia nitrogen concentration in the effluent, total nitrogen concentration in the effluent, total suspended solids concentration in the effluent, bio-oxygen demand concentration in the effluent, and chemical oxygen demand concentration in the effluent; the key process variable data were collected by sensors at time t as follows:
[0074] S(t)=[S O (t),S NO (t),S eNH (t),S eTN (t),S eTSS (t),S eBOD (t),S eCOD (t)] (1)
[0075] Where S(t) is the key variable vector at time t, S O (t) represents the dissolved oxygen concentration in the aerobic reactor at time t, in mg / L, S. NO (t) represents the nitrate nitrogen concentration in the anoxic reactor at time t, in mg / L, S eNH (t) represents the ammonia nitrogen concentration in the effluent at time t, in mg / L, S eTN (t) represents the total nitrogen concentration in the effluent at time t, in mg / L, S eTSS (t) represents the total suspended solids concentration in the effluent at time t, in mg / L, S eBOD (t) represents the biological oxygen demand at time t, in mg / L, S eCOD (t) represents the chemical oxygen demand at time t, in mg / L;
[0076] (2) Establish a time-delay robust multi-objective optimization model for urban wastewater treatment process: Use kernel functions to establish a time-delay robust multi-objective optimization model for urban wastewater treatment process, including effluent quality optimization model, operation energy consumption optimization model and time-delay robust compensation model, to express the relationship between key process variables and effluent quality and operation energy consumption;
[0077] ① Establish an effluent water quality optimization model:
[0078]
[0079] Among them, J EQ(x(t)) is the effluent water quality optimization model at time t, x(t) = [S O (t),S NO [t] represents J at time t. EQ The input vector of (x(t)), δ EQ,k (t) represents J at time t. EQ The k-th weight of (x(t)) and its value range is [-10, 10], c EQ,k (t)=[c EQ,k,1 (t),c EQ,k,2 [t] represents J at time t. EQ The k-th center vector of (x(t)), c EQ,k,1 (t) is c EQ,k The first component of (t) takes values in the range [-5, 5], c EQ,k,2 (t) is c EQ,k The second component of (t) takes values in the range [-5, 5], b EQ,k (t) represents J at time t. EQ The k-th width of (x(t)) and its value range is [0,3];
[0080] ② Establish an operational energy consumption optimization model:
[0081]
[0082] Among them, J EC (x(t)) is the energy consumption optimization model at time t, x(t) = [S O (t),S NO [t] represents J at time t. EC The input vector of (x(t)), δ EC,k (t) represents J at time t. EC The k-th weight of (x(t)) and its value range is [-10, 10], c EC,k (t)=[c EC,k,1 (t),c EC,k,2 [t] represents J at time t. EC The k-th center vector of (x(t)), c EC,k,1 (t) is c EC,k The first component of (t) takes values in the range [-5, 5], c EC,k,2 (t) is c EC,k The second component of (t) takes values in the range [-5, 5], b EC,k (t) represents J at time t. EC The k-th width of (x(t)) and its value range is [0,3];
[0083] ③ Establish a time-delay robust compensation model:
[0084]
[0085] Wherein, ε(x) τ (t) represents the time-delay robust compensation model at time t, x τ (t)=[S O (t-τ DO ),S NO (t-τ NO )] represents ε(x) at time t τ The input vector of τ(t)) DO =0.5(τ) DO,low +τ DO,up τ represents the time lag of dissolved oxygen in the aerobic reactor relative to the effluent. DO,low =V S / Q in For τ DO The lower bound, V S Q is the volume of the secondary sedimentation tank. in τ is the volume of water entering the vessel. DO,up =(V S +V O ) / Q in For τ DO The upper bound of V O S is the volume of the aerobic reactor. O (t-τ DO ) is t-τ DO The dissolved oxygen concentration in the aerobic reactor at any given time is within the range of [0,3], and the unit is mg / L; S NO (t-τ NO ) is t-τ NO The nitrate nitrogen concentration in the aerobic reactor at a given time is within the range of [0,2], and the unit is mg / L. NO =0.5(τ) DO,low +τ DO,up τ represents the time lag of nitrate nitrogen relative to the effluent in the anoxic reactor. NO,low =(V S +3V O ) / Q in For τ NO The lower bound of τ NO,up =(V S +V A +3V O ) / Q in For τ NO The upper bound of V A The volume of the anoxic reactor; δ ε,k (t) represents ε(x) at time t. τ The k-th weight of (t) is in the range [-10, 10], c ε,k (t)=[c ε,k,1(t),c ε,k,2 [(t)] represents ε(x) at time t. τ The k-th center vector of (t) is c ε,k,1 (t) is c ε,k The first component of (t) takes values in the range [-5, 5], c ε,k,2 (t) is c ε,k The second component of (t) takes values in the range [-5, 5], b ε,k (t) represents ε(x) at time t. τ The k-th width of (t) is in the range of [0,3].
[0086] ④ For J EQ (x(t)), J EC (x(t)) and ε(x) τ The parameters in (t) are updated in real time:
[0087]
[0088]
[0089]
[0090]
[0091]
[0092]
[0093]
[0094]
[0095]
[0096] Where, δ EQ (t)=[δ EQ,1 (t),δ EQ,2 (t),…,δ EQ,5 [t] represents J at time t. EQ The weight vector of (x(t)), δ EC (t)=[δ EC,1 (t),δ EC,2 (t),…,δ EC,5 [t] represents J at time t. EC The weight vector of (x(t)), δ ε (t)=[δ ε,1 (t),δ ε,2 (t),…,δ ε,5 [(t)] represents ε(x) at time t.τ The weight vector of (t)), c EQ (t)=[c EQ,1 (t),c EQ,2 (t),…,c EQ,5 [t] represents J at time t. EQ The center vector of (x(t)), c EC (t)=[c EC,1 (t),c EC,2 (t),…,c EC,5 [t] represents J at time t. EC The center vector of (x(t)), c ε (t)=[c ε,1 (t),c ε,2 (t),…,c ε,5 [(t)] represents ε(x) at time t. τ The center vector of (t)), b EQ (t)=[b EQ,1 (t),b EQ,2 (t),…,b EQ,5 [t] represents J at time t. EQ The width vector of (x(t)), b EC (t)=[b EC,1 (t),b EC,2 (t),…,b EC,5 [t] is J EC The width vector of (x(t)) at time t, b ε (t)=[b ε,1 (t),b ε,2 (t),…,b ε,5 [(t)] represents ε(x) at time t. τ The width vector of (t)), δ EQ (t+1)=[δ EQ,1 (t+1),δ EQ,2 (t+1),…,δ EQ,5 [t+1] represents J at time t+1. EQ The weight vector of (x(t+1)), δ EC (t+1)=[δ EC,1 (t+1),δ EC,2 (t+1),…,δ EC,5 [t+1] represents J at time t+1. EC The weight vector of (x(t+1)), δ ε (t+1)=[δ ε,1 (t+1),δ ε,2 (t+1),…,δ ε,5 [(t+1)] represents ε(x) at time t+1. τThe weight vector of (t+1)), c EQ (t+1)=[c EQ,1 (t+1),c EQ,2 (t+1),…,c EQ,5 [t+1] represents J at time t+1. EQ The center vector of (x(t+1)), c EC (t+1)=[c EC,1 (t+1),c EC,2 (t+1),…,c EC,5 [t+1] represents J at time t+1. EC The center vector of (x(t+1)), c ε (t+1)=[c ε,1 (t+1),c ε,2 (t+1),…,c ε,5 [(t+1)] represents ε(x) at time t+1. τ The center vector of (t+1)), b EQ (t+1)=[b EQ,1 (t+1),b EQ,2 (t+1),…,b EQ,5 [t+1] represents J at time t+1. EQ The width vector of (x(t+1)), b EC (t+1)=[b EC,1 (t+1),b EC,2 (t+1),…,b EC,5 [t+1] represents J at time t+1. EC The width vector of (x(t+1)), b ε (t+1)=[b ε,1 (t+1),b ε,2 (t+1),…,b ε,5 [(t+1)] represents ε(x) at time t+1. τ (t+1)) width vector; e EQ (t)=J EQ (x(t))-EQ(t) represents the training error of the effluent water quality optimization model at time t, EQ(t) = S eNH (t)+S eTN (t)+S eTSS (t)+S eBOD (t)+S eCOD (t) represents the effluent water quality sample value at time t, e EC (t)=J EC (x(t))-EC(t) represents the training error of the energy consumption optimization model at time t, and EC(t) represents the energy consumption collected at time t. ε (t)=e EQ(t)-ε(x τ (t) represents the time-delay robust compensation model training error at time t;
[0097] ⑤ Establish a time-delay robust multi-objective optimization function for the urban wastewater treatment process:
[0098]
[0099] J EQ (x(t)) represents the objective function value of the decision variable x(t) at time t, J EC (x(t)) represents the objective function value of the decision variable x(t) at time t, ε(x τ (t) represents the decision variable x. τ (t) The time delay robust compensation value at time t;
[0100] (3) Solving for optimal setpoints of process variables: Designing a time-delay robust evolutionary optimization algorithm to solve for optimal setpoints of dissolved oxygen and nitrate nitrogen;
[0101] ① Set the total number of iterations for the time-delay robust evolutionary optimization algorithm to κ. max =500, the particle swarm size is N=100, and the particle positions and velocities are:
[0102] z t,n (κ)=[z t,n,1 (κ),z t,n,2 (κ)](15)
[0103] v t,n (κ)=[v t,n,1 (κ),v t,n,2 (κ)](16) where z t,n (κ) is the position vector of the nth particle at time t, which evolves to the κth generation. t,n,1 (κ) is z t,n The first component of (κ) represents the dissolved oxygen concentration, and its value ranges from [0,3]. t,n,2 (κ) is z t,n The second component of (κ) represents the nitrate nitrogen concentration, with a value range of [0,2]. t,n (κ) is the velocity vector of the nth particle evolving to the κth generation at time t, v t,n,1 (κ) is v t,n The first component of (κ) and its value range is [0, 0.1], v t,n,2 (κ) is v t,n The second component of (κ) and its value range is [0, 0.1];
[0104] ② Calculate the individual optimal position of each particle and the global optimal position of the population:
[0105]
[0106]
[0107] Where, p t,n (κ+1) represents the optimal position of the nth particle at time t in the (κ+1)th generation, p t,n (κ) represents the optimal position of the nth particle at time t in the historical individual evolution to the κth generation, z t,n (κ+1) is the position vector of the nth particle evolving to the (κ+1)th generation at time t, and g t (κ+1) represents the globally optimal position of the population at time t after evolving to the (κ+1)th generation, and g t (κ) represents the globally optimal position of the population at time t in the κth generation, and > represents a robust dominance relationship defined as:
[0108]
[0109]
[0110] Among them, J EQ (z t,n (κ+1) is the decision variable z. t,n The objective function value for optimizing effluent water quality is (κ+1); J EC (z t,n (κ+1) is the decision variable z. t,n The objective function value for optimizing the operating energy consumption of (κ+1); ε(z) t,n (κ+1) is the decision variable z. t,n (κ+1) time-delay robust compensation value; J EQ (p t,n (κ) is the decision variable p t,n The objective function value for optimizing effluent water quality (κ); J EC (p t,n (κ) is the decision variable p t,n The objective function value for optimizing the operating energy consumption of (κ); ε(p t,n (κ)) is the decision variable (p) t,n (κ)) time delay robust compensation value; J EQ (g t,n (κ) is the decision variable g t,n The objective function value for optimizing effluent water quality (κ); J EC (g t,n (κ) is the decision variable g t,n The objective function value for optimizing the operating energy consumption of (κ); ε(g) t,n (κ)) is the decision variable (g) t,n(κ)) time delay robust compensation value; The expected value of the time-delay robust compensation model in the robust region:
[0111]
[0112] Ω=(S NO,up (t)-S NO,low (t))(S O,up (t)-S O,low (t)) (22)
[0113]
[0114]
[0115]
[0116]
[0117] Wherein, ε(z) t,n (κ) is the decision variable z t,n The time-delay robust compensation value of (κ), S NO,up (t) represents S at time t. NO The upper bound of (t), S NO,low (t) represents S at time t. NO The lower bound of S(t) O,up (t) represents S at time t. O The upper bound of (t), S O,low (t) represents S at time t. O (t) is the lower bound of , where Ω is the area of the robust region and λ is the sampling time;
[0118] ③ Update the particle's velocity and position:
[0119] z t,n (κ+1)=z t,n (κ)+v t,n (κ+1) (27)
[0120] v t,n (κ+1)=0.7v t,n (κ)+0.5μ1(p t,n (κ)-v t,n (κ))+0.5μ2(g t (κ)-v t,n (κ)) (28)
[0121] Among them, v t,n(κ+1) is the velocity vector of the nth particle evolving to the κ+1th generation at time t, μ1 is the random number of individual evolution with a value range of [0,1], and μ2 is the random number of population evolution with a value range of [0,1].
[0122] ④ Check if the evolutionary process has reached the stopping condition: If the evolutionary generation κ = 500, terminate the evolutionary process and jump to step ⑤; otherwise, increment the evolutionary generation κ by 1 and return to step ②.
[0123] ⑤ From the κth max The solution set Z of the generation t (κ max )=[z t,1 (κ max ),z t,2 (κ max ),…,z t,N (κ max A solution x is randomly selected from [ ] ] * (t)=[S * O (t),S * NO [(t)] serves as the optimized setpoint for dissolved oxygen and nitrate nitrogen, where S * O (t) represents the optimal setpoint for dissolved oxygen in the aerobic reactor at time t, S * NO (t) represents the optimal setpoint for nitrate nitrogen in the anoxic reactor at time t;
[0124] (4) Tracking control optimization setpoints:
[0125] ① Solve for the optimal setpoint x using a multivariable proportional-integral-derivative control algorithm. * (t) Perform tracking control:
[0126]
[0127] Where, Δu(t)=[ΔK L a(t), ΔQ a (t)] T Let ΔK be the matrix of operations at time t. L a(t) represents the change in the dissolved oxygen transfer coefficient at time t, ΔQ a (t) represents the change in circulating flow rate at time t, K p H is a proportionality constant with a value range of [0,5]. l H is the integral coefficient and its value ranges from [0,5]. d The differential coefficients are in the range [0,5], e(t) = x * (t)-x(t) is the control error vector at time t;
[0128] ② Adjust the dissolved oxygen transfer coefficient and internal reflux flow rate using Δu(t):
[0129] K L a(t+1)=K L a(t)+ΔK L a(t) (30)
[0130] Q a (t+1)=Q a (t)+ΔQ a (t) (31)
[0131] Among them, K L a(t+1) is the dissolved oxygen transfer coefficient at time t+1, K L a(t) is the dissolved oxygen transfer coefficient at time t, Q a (t+1) represents the internal return flow rate at time t+1, Q a (t) represents the internal return flow rate at time t;
[0132] (5) K is transmitted via the host computer L a(t) and Q a The value of (t) is transmitted to the programmable logic controller (PLC), and the digital signal is transmitted to the frequency converter via a digital-to-analog converter (DAC). The frequency converter adjusts the frequency of the oxygen supply pump and the internal return pump, thereby adjusting the dissolved oxygen concentration to S. * O (t), the nitrate nitrogen concentration was adjusted to S * NO (t).
Claims
1. A time-delay robust optimization control method for urban wastewater treatment processes, characterized in that, Includes the following steps: (1) Data collection of key process variables: Seven variables were selected as key process variables: dissolved oxygen concentration in the aerobic reactor, nitrate nitrogen concentration in the anoxic reactor, ammonia nitrogen concentration in the effluent, total nitrogen concentration in the effluent, total suspended solids concentration in the effluent, bio-oxygen demand concentration in the effluent, and chemical oxygen demand concentration in the effluent; the key process variable data were collected by sensors at time t as follows: (1); Where S(t) is the key variable vector at time t, S O (t) represents the dissolved oxygen concentration in the aerobic reactor at time t, in mg / L, S. NO (t) represents the nitrate nitrogen concentration in the anoxic reactor at time t, in mg / L, S eNH (t) represents the ammonia nitrogen concentration in the effluent at time t, in mg / L, S eTN (t) represents the total nitrogen concentration in the effluent at time t, in mg / L, S eTSS (t) represents the total suspended solids concentration in the effluent at time t, in mg / L, S eBOD (t) represents the biological oxygen demand at time t, in mg / L, S eCOD (t) represents the chemical oxygen demand at time t, in mg / L; (2) Establish a time-delay robust multi-objective optimization model for urban wastewater treatment process: Use kernel functions to establish a time-delay robust multi-objective optimization model for urban wastewater treatment process, including effluent quality optimization model, operation energy consumption optimization model and time-delay robust compensation model, to express the relationship between key process variables and effluent quality and operation energy consumption; ① Establish an effluent water quality optimization model: (2); Among them, J EQ (x(t)) represents the objective function value of the decision variable x(t) at time t, x(t) = [S O (t),S NO [t] represents J at time t. EQ The input vector of (x(t)), δ EQ,k (t) represents J at time t. EQ The k-th weight of (x(t)) takes values in the range [-10, 10], c EQ,k (t)=[c EQ,k,1 (t), c EQ,k,2 [t] represents J at time t. EQ The k-th center vector of (x(t)), c EQ,k,1 (t) is c EQ,k The first component of (t) takes values in the range [-5, 5], c EQ,k,2 (t) is c EQ,k The second component of (t) takes values in the range [-5, 5], b EQ,k (t) represents J at time t. EQ The k-th width of (x(t)) is in the range [0, 3]. ② Establish an operational energy consumption optimization model: (3); Among them, J EC (x(t)) represents the objective function value of the decision variable x(t) at time t, x(t) = [S O (t), S NO [t] represents J at time t. EC The input vector of (x(t)), δ EC,k (t) represents J at time t. EC The k-th weight of (x(t)) takes values in the range [-10, 10], c EC,k (t)=[c EC,k,1 (t), c EC,k,2 [t] represents J at time t. EC The k-th center vector of (x(t)), c EC,k,1 (t) is c EC,k The first component of (t) takes values in the range [-5, 5], c EC,k,2 (t) is c EC,k The second component of (t) takes values in the range [-5, 5], b EC,k (t) represents J at time t. EC The k-th width of (x(t)) is in the range [0, 3]. ③ Establish a time-delay robust compensation model: (4); Wherein, ε(x) τ (t) represents the decision variable x. τ (t) The time-delay robust compensation value at time t, x τ (t)=[S O (t-τ DO ), S NO (t-τ NO )] represents ε(x) at time t τ The input vector of τ(t)) DO =0.5(τ DO,low +τ DO,up τ represents the time lag of dissolved oxygen in the aerobic reactor relative to the effluent. DO,low =V S / Q in For τ DO The lower bound, V S Q is the volume of the secondary sedimentation tank. in τ is the volume of water entering the vessel. DO,up =(V S +V O ) / Q in For τ DO The upper bound of V O S is the volume of the aerobic reactor. O (t-τ DO ) is t-τ DO The dissolved oxygen concentration in the aerobic reactor at any given time is [0, 3], in mg / L; S NO (t-τ NO ) is t-τ NO The nitrate nitrogen concentration in the aerobic reactor at any given time is [0, 2], in mg / L. NO =0.5(τ DO,low +τ DO,up τ represents the time lag of nitrate nitrogen relative to the effluent in the anoxic reactor. NO,low =(V S +3V O ) / Q in For τ NO The lower bound, τ NO,up =(V S +V A +3V O ) / Q in For τ NO The upper bound, V A The volume of the anoxic reactor; δ ε,k (t) represents ε(x) at time t. τ The k-th weight of (t) is in the range [-10, 10], c ε,k (t)=[c ε,k,1 (t),c ε,k,2 [(t)] represents ε(x) at time t. τ The k-th center vector of (t) is c ε,k,1 (t) is c ε,k The first component of (t) takes values in the range [-5, 5], c ε,k,2 (t) is c ε,k The second component of (t) takes values in the range [-5, 5], b ε,k (t) represents ε(x) at time t. τ The width of (t) is the k-th width and its value range is [0, 3]. ④ For J EQ (x(t)), J EC (x(t)) and ε(x) τ The parameters in (t) are updated in real time: (5); (6); (7); (8); (9); (10); (11); (12); (13); Where, δ EQ (t)=[δ EQ,1 (t), δ EQ,2 (t), …, δ EQ,5 [t] represents J at time t. EQ The weight vector of (x(t)), δ EC (t)=[δ EC,1 (t), δ EC,2 (t), …, δ EC,5 [t] represents J at time t. EC The weight vector of (x(t)), δ ε (t)=[δ ε,1 (t), δ ε,2 (t), …, δ ε,5 [(t)] represents ε(x) at time t. τ The weight vector of (t)), c EQ (t)=[c EQ,1 (t), c EQ,2 (t), …, c EQ,5 [t] represents J at time t. EQ The center vector of (x(t)), c EC (t)=[c EC,1 (t), c EC,2 (t), …, c EC,5 [t] represents J at time t. EC The center vector of (x(t)), c ε (t)=[c ε,1 (t), c ε,2 (t), …, c ε,5 [(t)] represents ε(x) at time t. τ The center vector of (t)), b EQ (t)=[b EQ,1 (t), b EQ,2 (t), …, b EQ,5 [t] represents J at time t. EQ The width vector of (x(t)), b EC (t)=[b EC,1 (t), b EC,2 (t), …, b EC,5 [t] is J EC (x(t)) is the width vector at time t, b ε (t)=[b ε,1 (t), b ε,2 (t),…, b ε,5 [(t)] represents ε(x) at time t. τ The width vector of (t)), δ EQ (t+1)=[δ EQ,1 (t+1), δ EQ,2 (t+1), …, δ EQ,5 [t+1] represents J at time t+1. EQ The weight vector of (x(t+1)), δ EC (t+1)=[δ EC,1 (t+1), δ EC,2 (t+1), …, δ EC,5 [t+1] represents J at time t+1. EC The weight vector of (x(t+1)), δ ε (t+1)=[δ ε,1 (t+1), δ ε,2 (t+1), …, δ ε,5 [(t+1)] represents ε(x) at time t+1. τ The weight vector of (t+1)), c EQ (t+1)=[c EQ,1 (t+1), c EQ,2 (t+1), …, c EQ,5 [t+1] represents J at time t+1. EQ The center vector of (x(t+1)), c EC (t+1)=[c EC,1 (t+1), c EC,2 (t+1), …, c EC,5 [t+1] represents J at time t+1. EC The center vector of (x(t+1)), c ε (t+1)=[c ε,1 (t+1), c ε,2 (t+1), …, c ε,5 [(t+1)] represents ε(x) at time t+1. τ The center vector of (t+1)), b EQ (t+1)=[b EQ,1 (t+1), b EQ,2 (t+1), …, b EQ,5 [t+1] represents J at time t+1. EQ The width vector of (x(t+1)), b EC (t+1)=[b EC,1 (t+1), b EC,2 (t+1), …, b EC,5 [t+1] represents J at time t+1. EC The width vector of (x(t+1)), b ε (t+1)=[b ε,1 (t+1), b ε,2 (t+1), …, b ε,5 [(t+1)] represents ε(x) at time t+1. τ (t+1)) width vector; e EQ (t)= J EQ (x(t))-EQ(t) represents the training error of the effluent water quality optimization model at time t, EQ(t)=S eNH (t)+S eTN (t)+S eTSS (t)+S eBOD (t)+S eCOD (t) represents the effluent water quality sample value at time t, e EC (t)=J EC (x(t))-EC(t) represents the training error of the energy consumption optimization model at time t, and EC(t) represents the energy consumption collected at time t. ε (t)=e EQ (t)-ε(x τ (t) represents the time-delay robust compensation model training error at time t; ⑤ Establish a time-delay robust multi-objective optimization function for the urban wastewater treatment process: (14); J EQ (x(t)) represents the objective function value of the decision variable x(t) at time t, J EC (x(t)) represents the objective function value of the decision variable x(t) at time t, ε(x τ (t) represents the decision variable x. τ (t) The time delay robust compensation value at time t; (3) Solving for optimal setpoints of process variables: Designing a time-delay robust evolutionary optimization algorithm to solve for optimal setpoints of dissolved oxygen and nitrate nitrogen; ① Set the total number of iterations for the time-delay robust evolutionary optimization algorithm to κ. max =500, the particle swarm size is N=100, and the particle positions and velocities are: (15); (16); Where z t,n (κ) is the position vector of the nth particle at time t, which evolves to the κth generation. t,n,1 (κ) is z t,n The first component of (κ) represents the dissolved oxygen concentration, and its value ranges from [0, 3]. t,n,2 (κ) is z t,n The second component of (κ) represents the nitrate nitrogen concentration, with values ranging from [0, 2]. t,n (κ) is the velocity vector of the nth particle evolving to the κth generation at time t, v t,n,1 (κ) is v t,n The first component of (κ) and its value range is [0, 0.1], v t,n,2 (κ) is v t,n The second component of (κ) and its value range is [0, 0.1]; ② Calculate the individual optimal position of each particle and the global optimal position of the population: (17); (18); Where, p t,n (κ+1) represents the optimal position of the nth particle at time t in the (κ+1)th generation, p t,n (κ) represents the optimal position of the nth particle at time t in the historical individual evolution to the κth generation, z t,n (κ+1) is the position vector of the nth particle evolving to the (κ+1)th generation at time t, and g t (κ+1) represents the globally optimal position of the population at time t after evolving to the (κ+1)th generation, and g t (κ) represents the globally optimal position of the population at time t after evolving to the κth generation. Represents a robust dominance relation and is defined as follows: (19); (20); Among them, J EQ (z t,n (κ+1) is the decision variable z. t,n The objective function value for optimizing effluent water quality is (κ+1); J EC (z t,n (κ+1) is the decision variable z. t,n The objective function value for optimizing the operating energy consumption of (κ+1); ε(z) t,n (κ+1) is the decision variable z. t,n (κ+1) time-delay robust compensation value; J EQ (p t,n (κ) is the decision variable p t,n The objective function value for optimizing effluent water quality (κ); J EC (p t,n (κ) is the decision variable p t,n The objective function value for optimizing the operating energy consumption of (κ); ε(p t,n (κ)) is the decision variable (p) t,n (κ)) time delay robust compensation value; J EQ (g t,n (κ) is the decision variable g t,n The objective function value for optimizing effluent water quality (κ); J EC (g t,n (κ) is the decision variable g t,n The objective function value for optimizing the operating energy consumption of (κ); ε(g) t,n (κ)) is the decision variable (g) t,n The time-delay robust compensation value of (κ)); The expected value of the time-delay robust compensation model in the robust region: (21); (22); ; Wherein, ε(z) t,n (κ) is the decision variable z t,n The time-delay robust compensation value of (κ), S NO,up (t) represents S at time t. NO The upper bound of (t), S NO,low (t) represents S at time t. NO The lower bound of S(t) O,up (t) represents S at time t. O The upper bound of (t), S O,low (t) represents S at time t. O (t) is the lower bound of , Ω is the area of the robust region, and λ is the sampling time; ③ Update the particle's velocity and position: (27); (28); Among them, v t,n (κ+1) is the velocity vector of the nth particle evolving to the κ+1th generation at time t, μ1 is the random number of individual evolution with a value range of [0, 1], and μ2 is the random number of population evolution with a value range of [0, 1]. ④ Check if the evolutionary process has reached the stopping condition: If the evolutionary generation κ = 500, terminate the evolutionary process and jump to step ⑤; otherwise, increment the evolutionary generation κ by 1 and return to step ②. ⑤ From the κth max The solution set Z of the generation t (κ max )=[z t,1 (κ max ), z t,2 (κ max ), …, z t,N (κ max A solution x is randomly selected from [ ] ] * (t)=[S * O (t), S * NO [(t)] serves as the optimized setpoint for dissolved oxygen and nitrate nitrogen, where S * O (t) represents the optimal setpoint for dissolved oxygen in the aerobic reactor at time t, S * NO (t) represents the optimal setpoint for nitrate nitrogen in the anoxic reactor at time t; (4) Tracking control optimization setpoints: ① Solve for the optimal setpoint x using a multivariable proportional-integral-derivative control algorithm. * (t) Perform tracking control: (29) Where, ∆u(t)=[ΔK L a(t), ∆Q a (t)] T Let ΔK be the matrix of operations at time t. L a(t) represents the change in the dissolved oxygen transfer coefficient at time t, ∆Q a (t) represents the change in circulating flow rate at time t, K p H is a proportionality constant with a value range of [0,5]. l H is the integral coefficient and its value ranges from [0, 5]. d Let x be the differential coefficient and take values in the range [0, 5]. Then e(t) = x * (t)-x(t) is the control error vector at time t; ② Adjust the dissolved oxygen transfer coefficient and internal reflux flow rate using ∆u(t): (30); (31); Among them, K L a(t+1) is the dissolved oxygen transfer coefficient at time t+1, K L a(t) is the dissolved oxygen transfer coefficient at time t, Q a (t+1) represents the internal return flow rate at time t+1, Q a (t) represents the internal return flow rate at time t; (5) K is transmitted via the host computer L a(t) and Q a The value of (t) is transmitted to the programmable logic controller (PLC), and the digital signal is transmitted to the frequency converter via a digital-to-analog converter (DAC). The frequency converter adjusts the frequency of the oxygen supply pump and the internal return pump, thereby adjusting the dissolved oxygen concentration to S. * O (t), the nitrate nitrogen concentration was adjusted to S * NO (t).
Citation Information
Patent Citations
Urban sewage treatment process optimization control method based on dynamic multi-objective particle swarm optimization algorithm
CN110161995A
Municipal sewage treatment real-time optimization control method based on operation process information
CN111399455A