CN114417631A - Irrigation area water transmission and distribution system modeling method based on observation data - Google Patents
Irrigation area water transmission and distribution system modeling method based on observation data Download PDFInfo
- Publication number
- CN114417631A CN114417631A CN202210095758.2A CN202210095758A CN114417631A CN 114417631 A CN114417631 A CN 114417631A CN 202210095758 A CN202210095758 A CN 202210095758A CN 114417631 A CN114417631 A CN 114417631A
- Authority
- CN
- China
- Prior art keywords
- water
- area
- flow
- channel
- ith
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 title claims abstract description 265
- 238000003973 irrigation Methods 0.000 title claims abstract description 68
- 230000002262 irrigation Effects 0.000 title claims abstract description 66
- 238000000034 method Methods 0.000 title claims abstract description 60
- 238000009826 distribution Methods 0.000 title claims abstract description 41
- 230000005540 biological transmission Effects 0.000 title claims abstract description 30
- 238000005094 computer simulation Methods 0.000 title claims abstract description 10
- 230000008859 change Effects 0.000 claims abstract description 48
- 230000008569 process Effects 0.000 claims abstract description 18
- 238000005457 optimization Methods 0.000 claims abstract description 15
- 238000013461 design Methods 0.000 claims abstract description 11
- 238000010845 search algorithm Methods 0.000 claims abstract description 5
- 238000004364 calculation method Methods 0.000 claims description 17
- 239000011159 matrix material Substances 0.000 claims description 16
- 238000003860 storage Methods 0.000 claims description 15
- 230000033228 biological regulation Effects 0.000 claims description 12
- 238000012549 training Methods 0.000 claims description 9
- 238000012937 correction Methods 0.000 claims description 6
- 230000001105 regulatory effect Effects 0.000 claims description 6
- 239000013598 vector Substances 0.000 claims description 6
- 238000012544 monitoring process Methods 0.000 claims description 5
- 230000009471 action Effects 0.000 claims description 4
- 230000003068 static effect Effects 0.000 claims description 4
- 230000001133 acceleration Effects 0.000 claims description 3
- 238000004422 calculation algorithm Methods 0.000 claims description 3
- 230000001276 controlling effect Effects 0.000 claims description 3
- 238000011156 evaluation Methods 0.000 claims description 3
- 230000005484 gravity Effects 0.000 claims description 3
- 238000005259 measurement Methods 0.000 claims description 3
- 210000002569 neuron Anatomy 0.000 claims description 3
- 238000004886 process control Methods 0.000 claims description 3
- 238000001363 water suppression through gradient tailored excitation Methods 0.000 claims description 3
- 238000003062 neural network model Methods 0.000 claims description 2
- 238000012546 transfer Methods 0.000 claims 1
- 230000009466 transformation Effects 0.000 claims 1
- 238000007726 management method Methods 0.000 abstract description 7
- 238000002474 experimental method Methods 0.000 description 2
- 230000001902 propagating effect Effects 0.000 description 2
- 238000004088 simulation Methods 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 238000011144 upstream manufacturing Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000018109 developmental process Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000011218 segmentation Effects 0.000 description 1
- 239000004575 stone Substances 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2113/00—Details relating to the application field
- G06F2113/08—Fluids
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- Algebra (AREA)
- Computing Systems (AREA)
- Fluid Mechanics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Mathematical Physics (AREA)
- Pure & Applied Mathematics (AREA)
- Flow Control (AREA)
Abstract
The invention discloses an irrigation area water transmission and distribution system modeling method based on observation data, which relates to the technical field of water transmission and distribution system modeling and water resource management, and comprises the following steps: acquiring the length of a water return area of each channel pool and the initial numerical value of the power wave velocity of each channel section; correcting the flow transmission process of the channel pseudo-uniform flow area, and establishing a system model under initial parameters; carrying out online optimization on implicit parameters of the system model by using a heuristic search algorithm; and (3) establishing a model to predict the water depth change process of the pseudo-uniform flow area, and obtaining an irrigation area water transmission and distribution system model in a discrete state equation form optimized on line according to observation data, so that the accuracy and reliability of the model are guaranteed. The method can be practically applied to the prediction of the water level change process of the irrigation area water delivery and distribution canal system and the design of an automatic control system thereof, can improve the prediction precision of the irrigation area water level, provides a more accurate system model for the automatic control of irrigation and water delivery and distribution, and provides a support for effectively improving the water resource utilization efficiency and the management level of the irrigation area.
Description
Technical Field
The invention relates to the technical field of modeling of a water transmission and distribution system and water resource management, in particular to a method for modeling a water transmission and distribution system of an irrigation area based on observation data.
Background
In the world, 69% of the total annual water resource utilization is used for agricultural irrigation, and the proportion is higher in arid regions, for example, the agricultural water accounts for 84% in 2020 of Ningxia Hui autonomous region of China. The open channel water delivery is still the most main mode for water resource distribution in irrigation areas in China, the automatic control technology of the irrigation channel system can effectively solve the problems of slow information transmission, difficult optimization scheduling, inaccurate water quantity control and the like in the traditional manual control, avoids the occurrence of over-irrigation and under-irrigation to the maximum extent, and is an important trend for agricultural irrigation development.
The system modeling is an automatically controlled foundation stone, and the modeling accuracy directly influences the performance of the automatic control system. Because the saint-wien equation system for describing the open channel unsteady flow is a nonlinear partial differential equation system, no analytic solution exists mathematically, and a certain degree of linearization strategy is mainly adopted in actual control. An Integral Delay Model (IDM) proposed by Schuurmans (1995) is an approximate Model based on saint-wien equation linearization, the Model takes a stable water flow state under a design working condition as a special solution of the equation, a coefficient containing a variable in an original saint-wien equation is replaced by a corresponding numerical value of the special solution, the coefficient can be regarded as a constant in a certain range, and thus a linearized system Model is obtained.
However, the practical application of the model in irrigation areas has the disadvantages that the lateral water taking along the channel distribution is not considered, which results in that the uniform flow condition described in the IDM model does not exist, the water level along the channel continuously changes along with the change of the water taking flow, and the part above the channel water return area is actually a pseudo-uniform flow area; the second is the factor causing the water level change of the backwater area, except the flow change transmitted from the upstream to the downstream, the water storage quantity change of the pseudo-uniform flow area also exists, which is not considered in the IDM model; thirdly, the reliability of the parameter calculation method is insufficient, two main parameters of the IDM model are the area of a backwater area and the lag time, and two methods can be obtained, wherein one method is an estimation formula, the other method is a step experiment, the horizontal plane assumption based on the IDM model is often not in accordance with an actual system, the accuracy is limited, the latter method cannot be realized in field application due to overhigh cost, the experiment can be basically carried out only on a simulation model, the obtained parameters are deviated from the actual system, and the calculated parameters are basically only suitable for preset working conditions, and the parameter reliability is basically invalid when the operation state of the irrigation area is switched between high and low water levels or more unknown inflow or unknown loss exists in channels.
At present, a system model for full channel control cannot well describe the situation that water is taken from a large and medium irrigation district channel system along the way, and is difficult to deal with the time-varying characteristic in the channel operation process, so how to solve the problems and improve the prediction precision of the water level of an irrigation district control point, and further improve the water resource utilization efficiency and the management level of the irrigation district is a technical problem to be solved urgently by the technical personnel in the field.
Disclosure of Invention
In view of the above, the invention provides a method for modeling a water delivery and distribution system of an irrigation area based on observation data, which can improve the prediction precision of the water level of a control point of the irrigation area, provide a more accurate system model for automatic control of irrigation and water delivery and distribution, be practically applied to prediction of the water level change process of the irrigation area and automatic control design of the irrigation process, and provide support for effectively improving the water resource utilization efficiency and management level of the irrigation area.
In order to achieve the above purpose, the invention provides the following technical scheme:
an irrigation area water transmission and distribution system modeling method based on observation data comprises the following steps:
acquiring initial values of the length of a water return area of each channel pool and the power wave velocity of each channel section as initial parameters;
correcting the flow transmission process of the channel pseudo-uniform flow area, and establishing a system model under initial parameters;
carrying out online optimization on the length of a water return area and the dynamic wave speed of implicit parameters of the system model to obtain an optimized system model;
and establishing an irrigation area water transmission and distribution system model in a discrete state equation form optimized on line according to the observation data.
Optionally, the obtaining of the length of the water return area of each channel pool and the initial value of the power wave velocity of each channel section specifically includes the following steps:
the method comprises the following steps of obtaining a physical parameter set of a transmission and distribution system of an irrigation area, specifically:
the check gate has the functions of separating and controlling the water flow movement, and n regulating check gates are arranged in the ditch system to form n ditch pool units; the physical parameter set to be acquired includes: the pile number positions of all regulation control gates and water intake gates and the length L of each ditch pooliBottom slope DiRoughness CiThe shape parameters of each cross section of the channel, such as the trapezoidal section, the bottom width biCoefficient of slope MiEtc.; channel running flow Q under design conditionsiWater level target H before gate of regulating control gatesi；
Calculating the length of the water return area, and segmenting the ditch pool, specifically:
for the ith channel pool, calculating the channel running flow QsiLower normal water depth HniTrial calculation can be carried out by adopting a Barplowski formula, wherein the water depth h is given, the overflowing area A and the hydraulic radius R are determined according to the shape of the section, the corresponding flow Q is calculated, and the value of the water depth h is adjusted until the calculated flow Q and the channel running flow Q are obtainedsiEqual, the water depth h is the channel running flow QsiLower normal water depth HniThe flow rate Q is calculated as follows:
in the formula: a is the flow area, m2(ii) a R is hydraulic radius, which is equal to the ratio of the flow area to the wet circumference, m; ciThe roughness of the ith ditch pond; diIs the bottom slope of the ith canal pond; x is an exponential parameter in the talent-talent coefficient; under the given section shape, after the water depth h is determined, the flow area A and the hydraulic radius R can be correspondingly obtained;
length L of backwater zone of ith canal pondbiThe calculation formula of (a) is as follows:
after the length of the water return area is determined, the influence of the flow change of the water taking gate in the water return area on the water level is basically equivalent, and the water taking gate and the water level are combined into
The part of the non-backwater area is called pseudo-uniform flow area, N is selectediThe pile number position of each water taking gate is taken as a sectional point, and the pseudo-uniform flow area is divided into Ni+1 stage, for the j stage, merging all the water intake lock flows in the stage into
Calculating the dynamic wave velocity of each channel section, the dynamic wave velocity V of the jth channel section in the ith channel pooli,jThe calculation formula of (a) is as follows:
Vi,j＝ci,j+vi,j (4)；
wherein ,ci,jThe speed m/s of the dynamic wave of the jth channel section in the ith channel pool in the static water is obtained; v. ofi,jThe water flow speed of the jth canal section in the ith canal pond is m/s; g is the acceleration of gravity, m/s2；Ai,jIs the flow area, m, of the jth channel section in the ith channel pool2；Bi,jThe water surface width m of the jth canal section in the ith canal pond; ri,jThe hydraulic radius of the jth canal section in the ith canal pond is m; x is the same as formula (2); ci,jThe roughness of the jth channel section in the ith channel pool is shown; di,jIs the bottom slope of the jth channel section in the ith channel pool.
Optionally, the flow transmission process of the channel pseudo-uniform flow area is corrected specifically as follows:
the change of the water storage quantity can affect the change of the water level of the control point, so that the water storage quantity S of the jth channel section in the ith channel pooli,jIs equivalent to the water intake flow of the water intake gate in the canal section
The intake water flow correction value after considering the process is calculated as follows:
wherein ,ΔSi,jThe water storage capacity change value of the jth channel section in the ith channel pool is obtained; the delta t is a discretized time step length which is determined by the frequency of measurement of real-time monitoring equipment of an irrigation area and is generally 2-10 minutes;
is the water storage area m of the jth canal section in the ith canal basin2；Δhi,jFor the water depth change value of the jth canal section in the ith canal pond and for the irrigation area with a complete monitoring system, the water level in front of each gate has real-time observation data, namely delta hi,jCan be directly obtained from the observed data.
Optionally, the method for establishing a system model under the initial parameters, namely establishing a system model in a state space equation form by using the length of the water return region and the power wave velocity as the initial parameters, specifically comprises the following steps:
determining the process control object of the water transmission and distribution system as the water level H before the regulating and controlling gate at the downstream of the ith canal pitiDifference y from set water level targeti＝Hi-HsiThen, the water level change rule of the water return area is as follows:
wherein ,
water storage area m of the water return area2Obtaining according to the shape parameters of the cross section and the water depth;
is the input flow of the head, m3/s；τi,0Delay time s for transmitting the flow change of the canal head to the water return area; n is a radical ofiThe number of water intake gates;
the corrected total water intake quantity m of the jth canal section in the ith canal basin3/s；τi,jTransmitting the delay time s of the flow change of the jth channel section in the ith channel pool to the water return area;
is the total water intake flow of the water return area, m3/s；
Is the discharge quantity m of the tail of the canal3/s；
Lag time τi,rThe calculation formula of (a) is as follows:
wherein ,li,jThe length m of the jth channel section in the ith channel pool;
discretizing equation (9) for a continuous system yields:
combining variables of all the channels in the whole channel system into a column vector form, taking the water level deviation, the variation and the historical flow process as state variables x (k), taking the regulation action of a regulation and control gate as control variables u (k), taking the water gate flow as disturbance variables d (k), taking the water level deviation as output variables y (k), and establishing a system model in a state space equation form as follows:
x(k+1)＝Gx(k)+Hu(k)+Zd(k) (12)；
y(k)＝Cx(k) (13)；
wherein G is a state matrix, H is a control matrix, Z is a disturbance matrix, and C is an output matrix.
Optionally, the obtaining of the optimized system model specifically includes the following steps:
the system model is evaluated based on data such as water level before the gate, gate-passing flow and the like observed in real time in an irrigation area, and the method specifically comprises the following steps:
selecting a prediction time domain p of a system, taking historical data at the moment k as input, predicting the water level deviation at the future moment k +1, k +2, the page:
merging the system vectors of the prediction time domain to obtain:
the formula for predicting the water level using the formulas (12) and (13) is as follows:
wherein, the expression of the system matrix is as follows:
will predict the result
Comparing with the observation result Y (k), obtaining the root mean square error RMSE of the water level prediction as a performance index of the system model, and evaluating the system model;
a set of backwater zone length parameters L is givenbAnd a dynamic wave velocity V parameter, and obtaining a set of system model in a state space equation form and performance indexes by the method;
utilizing a heuristic search algorithm to carry out comparison on the length parameter L of the water return areabAnd optimizing the power wave velocity V parameter to obtain an optimized system model.
Optionally, for the length parameter L of the backwater areabThe specific method for optimizing the kinetic wave velocity V parameter is a differential evolution algorithm, and specifically comprises the following steps:
setting specific values of a variation factor F, a cross factor E, a population scale M and a maximum iteration number W;
randomly generating an initial population according to a certain rule, namely, parameter combinations of the lengths of M groups of backwater areas and the dynamic wave speed, and randomly carrying out variation and cross operation on each individual in the initial population;
when the performance index RMSE corresponding to the variant and crossed individuals is superior to the individual before the variant and crossed, replacing the original individual with the variant and crossed individual, otherwise, still retaining the original individual, and completing one iteration after the operation is performed on each individual;
and repeating the operation until the population individuals meet a certain threshold condition or reach the maximum iteration number W, stopping optimization, and outputting a parameter combination with optimal performance, wherein the parameter combination is the optimized system model according to the corresponding formula (12) and formula (13).
Optionally, the establishing of the irrigation area water transmission and distribution system model in the form of the discrete state equation optimized on line according to the observation data specifically includes the following steps:
establishing a prediction model of the water depth change of the pseudo-uniform flow area to correct the system model, which specifically comprises the following steps:
adopting a radial basis function neural network model RBF, taking observation data with a certain length in the near term as a training sample, inputting channel section water level variation and flow variation with a certain space-time range, setting neuron quantity N of a hidden layer, an expansion constant S of a radial basis function and a root mean square error target RMSE, and training to obtain an RBF prediction model, wherein the output quantity is the water level variation of a predicted channel section at the next moment;
correcting the water depth change of the pseudo-uniform flow area obtained based on D (k) obtained by a water demand plan and an RBF prediction model, and substituting an equation (16) to obtain a system model with prediction capability;
and (3) iteratively repeating the steps by the system based on the online updating of the observation data set, inputting the latest observation data into a system model for evaluation and training by optimizing the length of the backwater area and the dynamic wave speed and predicting the water depth, and establishing the irrigation area water transmission and distribution system model in a discrete state equation form according to the online optimization of the observation data.
Compared with the prior art, the invention discloses a method for modeling the water conveying and distributing system of the irrigation area based on the observation data, and the method has the following beneficial effects:
(1) the intelligent optimization algorithm adopted by the invention has the characteristic of rapid convergence, only second-level processing time is needed for online system identification based on a large number of observation data sequences of the irrigation area, the method can be practically applied to water level change process prediction of the irrigation area water transmission and distribution canal system and automatic control system design thereof, the prediction precision of the irrigation area water level can be improved, a more accurate system model is provided for automatic control of irrigation transmission and distribution, and support is provided for effectively improving the water resource utilization efficiency and management level of the irrigation area;
(2) the invention utilizes real-time observation data of the irrigation area to carry out online optimization on the parameters of the system model, and the system model obtained at each time interval has better fitting to the actual state of the current system, so the model can effectively adapt to the nonlinearity and time-varying property of the water delivery and distribution system;
(3) the invention solves the problems of excessive simplification, difficult accurate acquisition of parameters, limited applicable working condition range and the like of the current system model controlled by all channels, can better describe the situation that a large and medium irrigation area has more water intakes, and can cope with the time-varying phenomenon in the channel operation process.
Drawings
In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the embodiments or the prior art will be briefly described below, it is obvious that the drawings in the following description are only embodiments of the present invention, and for those skilled in the art, other drawings can be obtained according to the provided drawings without creative efforts.
FIG. 1 is an overall flow chart of a modeling method of a water delivery and distribution system of an irrigation area in the invention;
FIGS. 2(a) -2 (d) are diagrams illustrating the effect of the embodiment of the method of the present invention.
Detailed Description
The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.
The embodiment of the invention discloses an irrigation area water transmission and distribution system modeling method based on observation data, wherein the overall flow chart of the method is shown in figure 1, and the method comprises the following steps:
1) the method comprises the following steps of obtaining initial values of the length of a water return area of each channel pool and the power wave velocity of each channel section as initial parameters:
1-1) acquiring a physical parameter set of a transmission and distribution water system of an irrigation area, which specifically comprises the following steps:
in the channel of this embodiment, there are 1 water inlet floodgate, 3 regulation and control floodgates, and 3 regulation and control floodgates form 3 ditch pond units, there are 4 lateral water intaking floodgates in each ditch pond. Specific physical parameters are shown in table 1.
TABLE 1 physical parameters of irrigation canal systems
1-2) calculating the length of a backwater area, and segmenting a ditch pool;
for the ith channel pool, calculating the channel running flow QsiLower normal water depth HniTrial calculation can be carried out by adopting a Barplowski formula, wherein the water depth h is given, the overflowing area A and the hydraulic radius R are determined according to the shape of the section, the corresponding flow Q is calculated, and the value of the water depth h is adjusted until the calculated flow Q and the channel running flow Q are obtainedsiEqual, the water depth h is the channel running flow QsiLower normal water depth HniThe flow rate Q is calculated as follows:
in the formula: a is the flow area, m2(ii) a R is hydraulic radius, which is equal to the ratio of the flow area to the wet circumference, m; ciThe roughness of the ith ditch pond; diIs the bottom slope of the ith canal pond; x is an exponential parameter in the talent-talent coefficient; under the given section shape, after the water depth h is determined, the flow area A and the hydraulic radius R can be correspondingly obtained. And calculating to obtain the normal water depths of the three canal ponds of 1.65m, 1.87m and 1.60m in sequence under the design working condition.
Length L of backwater zone of ith canal pondbiThe calculation formula of (a) is as follows:
the lengths of the water return areas of the three canal ponds are 1250m, 650m and 1000m in sequence. The pile number positions of 4 water taking gates in each ditch pool are selected as segmentation points, and the ditch pool is divided into 5 segments.
1-3) calculating the dynamic wave velocity of each channel section;
the dynamic wave velocity is the propagation velocity of the fluctuation of the water flow in the static water, and the propagation velocity in the channel is added with the flow velocity of the water flow, because the water flow of the irrigation system is usually slow flow, the dynamic wave can simultaneously propagate upstream and downstream. Because the water level is controlled to be generally positioned at the most downstream of the ditch pool, the speed of the power wave propagating downstream can be considered, namely the speed V of the power wave of the jth ditch section in the ith ditch pooli,jThe calculation formula of (a) is as follows:
Vi,j＝ci,j+vi,j (4)；
wherein ,ci,jIs as followsThe speed m/s of the dynamic wave of the jth channel section in the i channel pools in the static water; v. ofi,jThe water flow speed of the jth canal section in the ith canal pond is m/s; g is the acceleration of gravity, m/s2；Ai,jIs the flow area, m, of the jth channel section in the ith channel pool2；Bi,jThe water surface width m of the jth canal section in the ith canal pond; ri,jThe hydraulic radius of the jth canal section in the ith canal pond is m; x is the same as formula (2); ci,jThe roughness of the jth channel section in the ith channel pool is shown; di,jIs the bottom slope of the jth channel section in the ith channel pool.
The corresponding section parameters can be calculated by the water depth, and the water depth can be obtained by flow trial calculation under the design working condition. The estimated values of the power wave velocity of 15 channel sections of 3 channel pools are obtained by calculation and are shown in table 2 in the unit of m/s.
TABLE 2 Power wave velocity estimation for each channel segment
2) Correcting the flow transmission process of the channel pseudo-uniform flow area, and establishing a system model under initial parameters, which comprises the following specific steps:
2-1) the change of the water storage capacity of each channel section is equivalent to the change of the water taking flow in the channel section;
the water storage quantity S of the jth canal section in the ith canal pondi,jThe change of the water level of the control point is influenced, and the equivalent of the change of the water level of the control point can be regarded as the water intake flow of the water intake gate in the canal section
The intake water flow correction value after considering the process is calculated as follows:
wherein ,ΔSi,jThe water storage capacity change value of the jth channel section in the ith channel pool is obtained; Δ t is a discretized time step, which is determined by the frequency of measurement performed by real-time monitoring equipment in an irrigation area, and is selected to be 2 minutes in the embodiment;
is the water storage area m of the jth canal section in the ith canal basin2Under the condition that the section shape is known, the water depth can be calculated; Δ hi,jFor the water depth change value of the jth canal section in the ith canal pond and for the irrigation area with a complete monitoring system, the water level in front of each gate has real-time observation data, namely delta hi,jCan be directly obtained from the observed data.
2-2) establishing a system model in a state space equation form by taking the length of a backwater area and the power wave velocity as initial parameters;
the water level H before the regulation and control gate at the downstream of the ith canal pitiDifference y from the set water level targeti＝Hi-HsiIf the water level is the target of the process control of the water delivery and distribution system, the water level change rule of the water return area is as follows:
wherein ,
water storage area m of the water return area2Obtained according to the profile shape parameters and the water depth, which are 33875m in this embodiment in sequence through calculation2、11531m2、16800m2；
Is the input flow of the head, m3/s；τi,0Delay time s for transmitting the flow change of the canal head to the water return area; n is a radical ofiThe number of water intake gates in the non-return water area;
the corrected total water intake quantity m of the jth canal section in the ith canal basin3/s；τi,jTransmitting the delay time s of the flow change of the jth channel section in the ith channel pool to the water return area;
is the total water intake flow of the water return area, m3/s；
Is the discharge quantity m of the tail of the canal3/s；
Lag time τi,rThe calculation formula of (a) is as follows:
wherein ,li,jThe length m of the jth channel section in the ith channel pool;
discretizing equation (9) for a continuous system yields:
the lag times for the three ponds obtained by calculation and sorting are shown in table 3 in units of s.
TABLE 3 estimated value of delay time for propagation of flow change of each gate to the backwater zone
Because the 3 rd and 4 th water taking gates are positioned in the water return area under the design condition, the delay time of the flow change propagating to the water return area is 0, however, the length of the water return area can be changed along with the operation of the system, and the corresponding steps are not suitable to be omitted. Combining variables of all channel ponds of the whole channel system into a column vector form, taking a water level deviation, a variable quantity and a historical flow process as state variables x (k), taking a regulation action of a regulation and control gate as control variables u (k), taking water gate flow as disturbance variables d (k), taking the water level deviation as output variables y (k), and establishing a system model in a state space equation form as follows:
x(k+1)＝Gx(k)+Hu(k)+Zd(k) (12)；
y(k)＝Cx(k) (13)；
wherein G is a state matrix, H is a control matrix, Z is a disturbance matrix, and C is an output matrix. The flow of the water inlet gate and the regulating control gate are sequentially marked as Q1(k)～Q4(k) And the 12 water intake gates are sequentially marked as q1(k)～q12(k) In this embodiment, values of each system variable and the coefficient matrix are as follows:
3) the method comprises the following steps of carrying out online optimization on the length of a water return area and the speed of a power wave, which are implicit parameters of a system model, and specifically comprising the following steps:
3-1) evaluating a system model according to data such as water level before a gate, gate passing flow and the like observed in real time in an irrigation area;
the prediction time domain p of the system is selected, in this embodiment, p is selected to be 5, historical data at the time k is used as input, and the water level deviation at the future time k +1, k + 2.
Merging the system vectors of the prediction time domain, namely:
the formula for predicting the water level using the formulas (12) and (13) is as follows:
wherein, the expression of the system matrix is as follows:
will predict the result
And comparing with the observation result Y (k) to obtain the root mean square error RMSE of the water level prediction as a performance index of the system model.
3-2) carrying out online optimization on the length of a water return area and the dynamic wave speed of hidden parameters of the system model by using a heuristic search algorithm;
a set of backwater zone length parameters L is givenbAnd dynamic wave velocity V parameters, and a set of system models and performance indexes in the form of state space equations are obtained by the above method, so that a heuristic search algorithm can be used to optimize two sets of parameters, taking a Differential Evolution (DE) as an example, in this embodiment:
setting a variation factor F to be 0.5, a cross factor E to be 0.7, a population scale M to be 200 and a maximum iteration number W to be 20; then, carrying out random variation on the initial parameters within a value range to generate an initial population, namely 200 groups of parameter combinations of the length of a backwater area and the dynamic wave velocity, and carrying out random variation and cross operation on each individual in the initial population; when the performance index RMSE corresponding to the individual after the variation and the intersection is superior to the individual before the variation and the intersection, the new individual is used for replacing the original individual, otherwise, the original individual is still reserved, and after the operation is carried out on each individual, one iteration is completed; and finally, repeating the operation until the population individual meets the threshold condition RMSE <0.0001m or the maximum iteration number W is 20, stopping optimization, and outputting a parameter combination with optimal performance, wherein the parameter combination is the corresponding formula (12) and formula (13), namely the optimized system model.
4) The method comprises the following steps of establishing an irrigation area water transmission and distribution system model in a discrete state space equation form optimized on line according to observation data, and specifically comprising the following steps:
4-1) establishing a prediction model of the water depth change of the pseudo-uniform flow area;
when the control point water level is predicted, in the formula (16), u (k) is a control action in a prediction time domain and is a known quantity for the control model, d (k) is a water intake gate flow change after correction in the prediction time domain, and the water intake gate flow change before correction is a known quantity for the control model, but the future correction quantity delta h in the formula (8)i,jThe model is unknown, so a prediction model of the water depth change of the pseudo-uniform flow area needs to be established to correct the model, and the error is reduced.
According to the water flow change rule of the pseudo-uniform flow area, the water depth change of a certain channel section is the result of the transmission of the water level and flow change of the local channel section and several nearby channel sections within a period of time, a non-linear model can be established to fit the process, taking a Radial Basis Function (RBF) as an example, taking the observation data with a certain length in the near term as a training sample, the input quantity is the water level change quantity and the flow change quantity of the channel section in a certain time-space range, the output quantity is the water level change quantity of the predicted channel section at the next moment, in the embodiment, the number N of neurons of an implied layer is set to be 20, the expansion constant S of a radial basis Function is 0.1, and the target RMSE is 10-6Then, the RBF prediction model can be obtained through training; according to D (k) obtained by a water demand plan, correcting by combining the water depth change of the pseudo-uniform flow area obtained by RBF prediction, and substituting the formula (16) to obtain the water flow area with prediction capabilityThe system model of (1).
4-2) establishing an online optimized water delivery and distribution system model;
with the online updating of the observation data set, the system iteratively repeats the steps, the parameter optimization of the step 3-2) and the water depth prediction of the step 4-1), and the latest observation data is included in the model evaluation and training, so that the continuous online optimization of the system model with the updating of the real-time observation data is realized.
In the embodiment, on-line modeling is performed on 16 gates in 3 canal ponds, the dimension of a system state variable is 17, the dimension of a control variable is 3, the dimension of a disturbance variable is 20, a prediction time domain is 5, the modeling is realized by MATLAB programming, the test finds that only 1-2 seconds are needed for completing one-time modeling, the average time for each modeling is 1.37 seconds for 200 times of test, the requirement of real-time control on the modeling speed is completely met, and the method can be applied to actual production.
The comparison between the water level variation process predicted by the linearized water transportation and distribution system model established by the present invention and the nonlinear numerical simulation calculation result is shown in fig. 2(a) -2 (d), in this embodiment, the water intake flow of 12 water intake gates is 0.3m at the beginning3And/s, the water level of each channel pool control point is equal to a design value, the system is stable and unchanged, and the error of the system model is also kept to be 0. At 1:42, an extreme disturbance situation that a water intake gate is suddenly and completely closed occurs, the water depth of a control point suddenly rises, the prediction of the water level by the system model begins to have an error, and the prediction error basically does not exceed 1cm in the whole tracking process. From the root mean square error index, when the water level changes suddenly, the prediction error is the largest and is about 1.5cm, then the prediction error rapidly decreases, the prediction error is basically maintained at 0.1-0.2 cm under the condition that the water level changes continuously, the model precision is ideal, the method can be practically applied to prediction of the water level change process of the irrigation area and automatic control design of a gate, the prediction precision of the water level of the irrigation area is improved, a more accurate system model is provided for automatic control of irrigation, transportation and water distribution, and support is provided for effectively improving the water resource utilization efficiency and the management level of the irrigation area.
The previous description of the disclosed embodiments is provided to enable any person skilled in the art to make or use the present invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the generic principles defined herein may be applied to other embodiments without departing from the spirit or scope of the invention. Thus, the present invention is not intended to be limited to the embodiments shown herein but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.
Claims (7)
1. An irrigation area water transmission and distribution system modeling method based on observation data is characterized by comprising the following steps:
acquiring initial values of the length of a water return area of each channel pool and the power wave velocity of each channel section as initial parameters;
correcting the flow transmission process of the channel pseudo-uniform flow area, and establishing a system model under initial parameters;
carrying out online optimization on the length of a water return area and the dynamic wave speed of implicit parameters of the system model to obtain an optimized system model;
and establishing an irrigation area water transmission and distribution system model in a discrete state equation form optimized on line according to the observation data.
2. The method for modeling the irrigation area water delivery and distribution system based on the observation data as claimed in claim 1, wherein the step of obtaining the initial values of the length of the water return area of each channel and the power wave velocity of each channel section comprises the following steps:
the method comprises the following steps of obtaining a physical parameter set of a transmission and distribution system of an irrigation area, specifically:
forming n trench pool units according to n regulation and control gates in the trench system; the obtained physical parameter set comprises: the pile number positions of all regulation control gates and water intake gates and the length L of each ditch pooliBottom slope DiRoughness CiShape parameters of each cross section of the channel and channel running flow Q under design working conditionsiWater level target H before gate of regulating control gatesi；
Calculating the length of the water return area, and segmenting the ditch pool, specifically:
for the ith channel pool, the water depth h is given, the flow area A and the hydraulic radius R are determined according to the section shape, the corresponding flow Q is calculated, and the value of the water depth h is adjusted until the calculated flow Q and the channel running flow Q are obtainedsiEqual, the water depth h is the channel running flow QsiLower normal water depth HniThe flow rate Q is calculated as follows:
in the formula: a is the flow area, m2(ii) a R is hydraulic radius, which is equal to the ratio of the flow area to the wet circumference, m; ciThe roughness of the ith ditch pond; diIs the bottom slope of the ith canal pond; x is an exponential parameter in the talent-talent coefficient;
length L of backwater zone of ith canal pondbiThe calculation formula of (a) is as follows:
calculating the dynamic wave velocity of each channel section, the dynamic wave velocity V of the jth channel section in the ith channel pooli,jThe calculation formula of (a) is as follows:
Vi,j＝ci,j+vi,j (4)；
wherein ,ci,jthe speed m/s of the dynamic wave of the jth channel section in the ith channel pool in the static water is obtained; v. ofi,jThe water flow speed of the jth canal section in the ith canal pond is m/s; g is the acceleration of gravity, m/s2；Ai,jIs the flow area, m, of the jth channel section in the ith channel pool2；Bi,jThe water surface width m of the jth canal section in the ith canal pond; ri,jThe hydraulic radius of the jth canal section in the ith canal pond is m; x is the same as formula (2); ci,jThe roughness of the jth channel section in the ith channel pool is shown; di,jIs the bottom slope of the jth channel section in the ith channel pool.
3. The modeling method for the irrigation area water delivery and distribution system based on the observation data as claimed in claim 1, wherein the flow transfer process of the channel pseudo-uniform flow area is corrected by:
the water storage quantity S of the jth channel section in the ith channel pondi,jIs equivalent to the water intake flow of the water intake gate in the canal section
The intake water flow correction value after considering the process is calculated as follows:
wherein ,ΔSi,jThe water storage capacity change value of the jth channel section in the ith channel pool is obtained; the delta t is a discretized time step length and is determined by the frequency of measurement of real-time monitoring equipment of an irrigation area;
is the water storage area m of the jth canal section in the ith canal basin2；Δhi,jFor the jth canal section in the ith canal pitAnd the water depth change value is obtained from the observation data.
4. The irrigation area water delivery and distribution system modeling method based on observation data as claimed in claim 2, wherein the system model under initial parameters is established, that is, the system model in the form of a state space equation is established by taking the length of the backwater area and the dynamic wave speed as initial parameters, and the method specifically comprises the following steps:
determining the process control object of the water transmission and distribution system as the water level H before the regulating and controlling gate at the downstream of the ith canal pitiDifference y from set water level targeti＝Hi-HsiThen, the water level change rule of the water return area is as follows:
wherein ,
water storage area m of the water return area2Obtaining according to the shape parameters of the cross section and the water depth;
is the input flow of the head, m3/s；τi,0Delay time s for transmitting the flow change of the canal head to the water return area; n is a radical ofiThe number of water intake gates;
the corrected total water intake quantity m of the jth canal section in the ith canal basin3/s；τi,jTransmitting the delay time s of the flow change of the jth channel section in the ith channel pool to the water return area;
is the total water intake flow of the water return area, m3/s；
Is the discharge quantity m of the tail of the canal3/s；
Lag time τi,rThe calculation formula of (a) is as follows:
wherein ,li,jThe length m of the jth channel section in the ith channel pool;
discretizing formula (9) yields:
combining variables of all the channels in the whole channel system into a column vector form, taking the water level deviation, the variation and the historical flow process as state variables x (k), taking the regulation action of a regulation and control gate as control variables u (k), taking the water gate flow as disturbance variables d (k), taking the water level deviation as output variables y (k), and establishing a system model in a state space equation form as follows:
x(k+1)＝Gx(k)+Hu(k)+Zd(k) (12)；
y(k)＝Cx(k) (13)；
wherein G is a state matrix, H is a control matrix, Z is a disturbance matrix, and C is an output matrix.
5. The irrigation area water delivery and distribution system modeling method based on observation data as claimed in claim 4, wherein the obtaining of the optimized system model specifically comprises the following steps:
evaluating the system model based on the water level before the gate and the gate passing flow observed in real time in the irrigation area, which specifically comprises the following steps:
selecting a prediction time domain p of a system, taking historical data at the moment k as input, predicting the water level deviation at the future moment k +1, k +2, the page:
merging the system vectors of the prediction time domain to obtain:
the formula for predicting the water level using the formulas (12) and (13) is as follows:
wherein, the expression of the system matrix is as follows:
will predict the result
Comparing with the observation result Y (k), obtaining the root mean square error RMSE of the water level prediction as a performance index of the system model, and evaluating the system model;
a set of backwater zone length parameters L is givenbAnd a dynamic wave velocity V parameter, and obtaining a set of system model in a state space equation form and performance indexes by the method;
utilizing a heuristic search algorithm to carry out comparison on the length parameter L of the water return areabAnd optimizing the power wave velocity V parameter to obtain an optimized system model.
6. The method for modeling the water conveying and distributing system of the irrigation area based on the observation data as claimed in claim 5, wherein the length parameter L of the backwater area is determined by the following formulabOptimizing with the V parameter of the velocity of the power waveThe specific method of the transformation is a differential evolution algorithm, and specifically comprises the following steps:
setting specific values of a variation factor F, a cross factor E, a population scale M and a maximum iteration number W;
randomly generating an initial population according to a certain rule, namely, parameter combinations of the lengths of M groups of backwater areas and the dynamic wave speed, and randomly carrying out variation and cross operation on each individual in the initial population;
when the performance index RMSE corresponding to the variant and crossed individuals is superior to the individual before the variant and crossed, replacing the original individual with the variant and crossed individual, otherwise, still retaining the original individual, and completing one iteration after the operation is performed on each individual;
and repeating the operation until the population individuals meet a certain threshold condition or reach the maximum iteration number W, stopping optimization, and outputting a parameter combination with optimal performance, wherein the parameter combination is the optimized system model according to the corresponding formula (12) and formula (13).
7. The method for modeling the irrigation area water delivery and distribution system based on the observation data as claimed in claim 5, wherein the method for establishing the irrigation area water delivery and distribution system model in the form of the discrete state equation optimized on line according to the observation data comprises the following steps:
establishing a prediction model of the water depth change of the pseudo-uniform flow area to correct the system model, which specifically comprises the following steps:
adopting a radial basis function neural network model RBF, taking observation data with a certain length in the near term as a training sample, inputting channel section water level variation and flow variation with a certain space-time range, setting neuron quantity N of a hidden layer, an expansion constant S of a radial basis function and a root mean square error target RMSE, and training to obtain an RBF prediction model, wherein the output quantity is the water level variation of a predicted channel section at the next moment;
correcting the water depth change of the pseudo-uniform flow area obtained based on D (k) obtained by a water demand plan and an RBF prediction model, and substituting an equation (16) to obtain a system model with prediction capability;
and (3) iteratively repeating the steps by the system based on the online updating of the observation data set, inputting the latest observation data into a system model for evaluation and training by optimizing the length of the backwater area and the dynamic wave speed and predicting the water depth, and establishing the irrigation area water transmission and distribution system model in a discrete state equation form according to the online optimization of the observation data.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210095758.2A CN114417631B (en) | 2022-01-26 | 2022-01-26 | Modeling method of infusion area water delivery and distribution system based on observation data |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210095758.2A CN114417631B (en) | 2022-01-26 | 2022-01-26 | Modeling method of infusion area water delivery and distribution system based on observation data |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114417631A true CN114417631A (en) | 2022-04-29 |
| CN114417631B CN114417631B (en) | 2023-10-03 |
Family
ID=81276558
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210095758.2A Active CN114417631B (en) | 2022-01-26 | 2022-01-26 | Modeling method of infusion area water delivery and distribution system based on observation data |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114417631B (en) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115456422A (en) * | 2022-09-16 | 2022-12-09 | 中国水利水电科学研究院 | Irrigation district water distribution plan dynamic preview correction method based on computational hydrodynamics |
| CN116822180A (en) * | 2023-06-19 | 2023-09-29 | 中国水利水电科学研究院 | Integral time-lag model parameter identification method for simulating river and canal water delivery process |
| CN117151187A (en) * | 2023-08-17 | 2023-12-01 | 清华大学 | Gate opening prediction model training application method, device, equipment and storage medium |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2015094122A (en) * | 2013-11-12 | 2015-05-18 | 沖電気工業株式会社 | Water level prediction device, water level prediction method and water level prediction system |
| CN106874622A (en) * | 2017-03-13 | 2017-06-20 | 中国水利水电科学研究院 | One kind series connection Duo Qu ponds gate lock flow coefficient rating method |
| CN109197539A (en) * | 2018-10-25 | 2019-01-15 | 清华大学 | A kind of irrigation multistage canal automatic control method based on Model Predictive Control Algorithm |
| CN113361214A (en) * | 2021-06-17 | 2021-09-07 | 武汉大学 | Open channel control model parameter identification method based on water level flow data |
-
2022
- 2022-01-26 CN CN202210095758.2A patent/CN114417631B/en active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2015094122A (en) * | 2013-11-12 | 2015-05-18 | 沖電気工業株式会社 | Water level prediction device, water level prediction method and water level prediction system |
| CN106874622A (en) * | 2017-03-13 | 2017-06-20 | 中国水利水电科学研究院 | One kind series connection Duo Qu ponds gate lock flow coefficient rating method |
| CN109197539A (en) * | 2018-10-25 | 2019-01-15 | 清华大学 | A kind of irrigation multistage canal automatic control method based on Model Predictive Control Algorithm |
| CN113361214A (en) * | 2021-06-17 | 2021-09-07 | 武汉大学 | Open channel control model parameter identification method based on water level flow data |
Non-Patent Citations (2)
| Title |
|---|
| 王忠静 等: "基于线性二次型的多级联输水渠道最优控制", 《水科学进展》 * |
| 苏海旺 等: "基于ID预测模型的明渠系统控制参数在线优化研究", 《中国农村水利水电》 * |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115456422A (en) * | 2022-09-16 | 2022-12-09 | 中国水利水电科学研究院 | Irrigation district water distribution plan dynamic preview correction method based on computational hydrodynamics |
| CN115456422B (en) * | 2022-09-16 | 2024-02-02 | 中国水利水电科学研究院 | Dynamic previewing correction method for irrigation area water distribution plan based on computational hydrodynamics |
| CN116822180A (en) * | 2023-06-19 | 2023-09-29 | 中国水利水电科学研究院 | Integral time-lag model parameter identification method for simulating river and canal water delivery process |
| CN116822180B (en) * | 2023-06-19 | 2024-01-19 | 中国水利水电科学研究院 | Integral time-lag model parameter identification method for simulating river and canal water delivery process |
| CN117151187A (en) * | 2023-08-17 | 2023-12-01 | 清华大学 | Gate opening prediction model training application method, device, equipment and storage medium |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114417631B (en) | 2023-10-03 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN114417631A (en) | Irrigation area water transmission and distribution system modeling method based on observation data | |
| CN110984062A (en) | Simulation scheduling method for large-scale reservoir group of watershed main and branch flows | |
| CN103730006B (en) | A kind of combination forecasting method of Short-Term Traffic Flow | |
| Unver et al. | Model for real-time optimal flood control operation of a reservoir system | |
| CN109711617B (en) | Medium-and-long-term runoff prediction method based on BLSTM deep learning | |
| CN111474965B (en) | Fuzzy neural network-based method for predicting and controlling water level of series water delivery channel | |
| Li et al. | Identifying explicit formulation of operating rules for multi-reservoir systems using genetic programming | |
| Bayat et al. | Optimization–simulation for short-term reservoir operation under flooding conditions | |
| CN108985577B (en) | Reservoir group real-time flood control dispatching significant reservoir intelligent identification method based on inference machine | |
| Ferdowsi et al. | Hybrid bat & particle swarm algorithm for optimization of labyrinth spillway based on half & quarter round crest shapes | |
| Salmasi et al. | Neuro-fuzzy approach for estimating energy dissipation in skimming flow over stepped spillways | |
| CN112799310A (en) | Method for urban drainage system simulation control mixed model based on mechanism model, concept model and data model of C language | |
| CN116523269B (en) | Model predictive control-based river network dynamic intelligent regulation and control method | |
| CN112596386A (en) | Matlab-based urban drainage system simulation control mixed model with mechanism model, concept model and data model | |
| CN112464584A (en) | Method for estimating water level and flow of free surface flow | |
| Fatih | Prediction of dam reservoir volume fluctuations using adaptive neuro fuzzy approach | |
| Kong et al. | Water level predictive control in multi-cascade pools with large control interval demand based on constraint methods | |
| Shahverdi et al. | Double Q-PI architecture for smart model-free control of canals | |
| CN113361214A (en) | Open channel control model parameter identification method based on water level flow data | |
| CN113050430A (en) | Drainage system control method based on robust reinforcement learning | |
| Yan et al. | A novel IBAS-ELM model for prediction of water levels in front of pumping stations | |
| CN101899820A (en) | Method for determining amount of available surface water of river basin facing to protection of river ecosystem | |
| Anilan et al. | Application of Artificial Neural Networks and regression analysis to L-moments based regional frequency analysis in the Eastern Black Sea Basin, Turkey | |
| Mu et al. | Real-time flood forecasting method with 1-D unsteady flow model | |
| Mohammad Rezapour Tabari et al. | Development of GWO–DSO and PSO–DSO hybrid models to redesign the optimal dimensions of labyrinth spillway |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |