A parameter identification method for improving the polarization characteristic model of proton exchange membrane electrolyzer

Abstract

The present invention discloses a parameter identification method for improving the polarization characteristic model of a proton exchange membrane electrolyzer, which belongs to the field of system modeling and parameter identification, and the steps are: establishing a polarization characteristic model of a proton exchange membrane electrolyzer and determining the parameters of the model to be identified; obtaining a measurable operating variable data set of the proton exchange membrane electrolyzer, including a measurable input and output variable data set; based on an adaptive hunting game strategy, performing parameter identification of the polarization characteristic model of the proton exchange membrane electrolyzer, and then obtaining the final polarization characteristic model of the proton exchange membrane electrolyzer. The present invention improves the polarization characteristic model in combination with the operating mechanism of the electrolyzer, and improves the modeling accuracy; proposes a parameter identification method based on an adaptive hunting game strategy, which, under the premise of ensuring accurate optimization of parameters, achieves rapid convergence of the parameter identification process and reduces the number of iterations required for parameter identification under a small population.

Description

Technical Field

[0001] The invention belongs to the technical field of proton exchange membrane electrolyzer modeling and parameter identification, and specifically relates to a parameter identification method for improving a polarization characteristic model of a proton exchange membrane electrolyzer. Background Art

[0002] With the continuous development of social economy, energy demand has also increased rapidly. The extensive use of fossil energy has aggravated climate change and environmental pollution. In this context, hydrogen energy is widely regarded as the most promising secondary energy in the future due to its green and pollution-free nature. As an efficient and environmentally friendly water electrolysis hydrogen production equipment, proton exchange membrane electrolyzer is increasingly widely used in hydrogen energy production. An accurate polarization characteristic model of proton exchange membrane electrolyzer is an important basis for achieving observable and controllable key states of proton exchange membrane electrolyzer, which requires accurate modeling and parameter identification of the polarization characteristics of the electrolyzer.

[0003] At present, a model of the polarization characteristics of proton exchange membrane electrolyzers has been proposed, but this model is inherited from the research on proton exchange membrane fuel cells and has not been modified for the operating mechanism of proton exchange membrane electrolyzers. It is often unable to accurately describe the actual working state of the electrolyzer. For the identification of unknown parameters in the polarization characteristics model of proton exchange membrane electrolyzers, existing naturally inspired swarm intelligence methods, such as particle swarm optimization methods and ant colony optimization methods, are prone to fall into local optimality, have slow convergence speed, and have poor parameter identification effects under small population sizes. They are not suitable for fast and accurate parameter identification of industrial systems with limited performance and resource allocation. Summary of the invention

[0004] In view of the above technical problems existing in the prior art, the present invention provides a parameter identification method for improving the polarization characteristic model of a proton exchange membrane electrolyzer, which has a reasonable design, solves the shortcomings of the prior art and has good effects.

[0005] In order to achieve the above object, the present invention adopts the following technical scheme:

[0006] A parameter identification method for improving a polarization characteristic model of a proton exchange membrane electrolyzer comprises the following steps:

[0007] S1. Establish a polarization characteristic model of a proton exchange membrane electrolyzer and determine the model parameters to be identified;

[0008] S2. Obtain a measurable operating variable data set of the proton exchange membrane electrolyzer, including a measurable input and output variable data set;

[0009] S3. Based on the adaptive hunting game strategy, the parameters of the polarization characteristic model of the proton exchange membrane electrolyzer are identified to obtain the final polarization characteristic model of the proton exchange membrane electrolyzer.

[0010] Furthermore, in step S1, the established proton exchange membrane electrolyzer polarization characteristic model is specifically:

[0011] The proton exchange membrane electrolyzer is composed of several units, and the working voltage of each unit is u cell Expressed as open circuit voltage u ocv , activation overvoltage η act 、Ohmic overvoltage η ohm Concentration difference overvoltage η conc The sum of the four is shown in the following formula:

[0012] u cell =u ocv +η act +η ohm +η conc (1)

[0013] Open circuit voltage u ocv The expression is:

[0014]

[0015] Where T is the thermodynamic temperature scale, the unit is K, measured by experiment; R is the gas constant; F is the Faraday constant; k wo is the correction factor; p H2 is the hydrogen pressure in the electrolyzer, p O2 is the oxygen pressure in the electrolytic cell, and its expression is:

[0016] p H2 =P ca -p H2O (3)

[0017] p O2 =P an -p H2O (4)

[0018] Among them, P ca and P an are the absolute pressures measured at the hydrogen separator and the oxygen separator, respectively, p H2O is the water pressure in the electrolytic cell, and its expression is:

[0019]

[0020] Activation overvoltage η act The expression is:

[0021]

[0022] Among them, α an is the anode charge transfer coefficient, i is the current flowing through each electrolytic unit, i 0,anis the anode current density, and its expression is:

[0023]

[0024] Among them, i 0,an,std is the anode current density under standard conditions, E exc The activation energy required for electron transport to the anode;

[0025] Ohmic overvoltage η ohm The expression is:

[0026]

[0027] Among them, σ mem is the conductivity of the proton exchange membrane, and its expression is:

[0028]

[0029] Among them, σ mem,std is the conductivity of the proton exchange membrane under standard conditions, E pro The activation energy required for proton transport in the membrane;

[0030] Since the rated operating current density in the experiment is not high enough to produce concentration difference overvoltage η conc , so the value is 0.

[0031] The six model parameters that need to be identified in the polarization characteristic model of the proton exchange membrane electrolyzer are: correction coefficient k wo , anode charge transfer coefficient α an , the anode current density i under standard conditions 0,an,std , the activation energy E required for anode electron transfer exc , membrane resistivity σ under standard conditions mem,std and the activation energy E required for proton transport in the membrane pro .

[0032] Furthermore, in step S2, experimental data of measurable parameters of the electrolyzer and corresponding working status of the electrolyzer are obtained and recorded through actual working condition testing of the proton exchange membrane electrolyzer, and the measurable parameters include voltage, current, temperature and pressure when the electrolyzer is operating.

[0033] Furthermore, the step S3 includes the following sub-steps:

[0034] S301, determine the reasonable range of the model parameters to be identified, and randomly generate initial population individuals within the range, specifically:

[0035] The range of the model parameters to be identified is given by existing public data or based on experience, and its expression is:

[0036]

[0037] Among them, lb1, lb2, lb3, lb4, lb5 and lb6 represent the lower bounds of the parameters to be identified, and ub1, ub2, ub3, ub4, ub5 and ub6 represent the upper bounds of the parameters to be identified;

[0038] According to the reasonable range of the model parameters to be identified, the initial population individuals are randomly generated, and the population number is set to N, that is, the number of individuals in the population. Each individual is a candidate solution, containing variables of 6 dimensions, representing the 6 parameters that need to be identified; the variable expression of the dth dimension of each candidate solution is:

[0039] x d =lb d +(ub d -lb d )·rand,d=1,2,...,6 (11)

[0040] Among them, x d represents the variable of the dth dimension, rand represents a random number in the range of [0,1], which is randomly determined when generating the variables of each dimension of all individuals;

[0041] S302. Evaluate individual danger intensity and competition intensity, specifically:

[0042] The root mean square error RMSE between the measured output and the model output is used as the fitness function f fitness (x k ), calculate the fitness of each individual, the expression is:

[0043]

[0044]

[0045] Among them, x k represents the position of the kth individual, x k =[x 1 ,x 2 ,x 3 ,x 4 ,x 5 ,x 6 ], h is the number of sample points contained in the measured data sequence, y i is the measured voltage output of the i-th sample point, is the model voltage output corresponding to the i-th sample point;

[0046] The individual with the lowest fitness is regarded as prey, and the other individuals in the population are regarded as predators;

[0047] For prey, DIprey It is defined as the intensity of danger perceived by the prey individual, and its expression is:

[0048]

[0049] Where m represents the number of predators that the prey can perceive, which is determined by roulette; x prey and x hunter,m Represent the positions of the prey and the mth predator at the perceptible boundary respectively; r1 and r2 are two random numbers in the range [0,1], and D represents the number of dimensions of the candidate solution, that is, the number of parameters. is the variable of the d-th dimension of the prey, is the variable of the jth predator in the dth dimension that the prey can perceive, α represents the endurance performance of the prey, which decreases with the increase of the number of iterations, and the expression is:

[0050]

[0051] Among them, C is the endurance constant, t is the current number of iterations, and t max The maximum number of iterations set;

[0052] For predators, DI hunter It is defined as the competition intensity perceived by individual predators and is expressed as:

[0053]

[0054] Where n represents the number of surrounding predators that the predator can perceive, which is determined by roulette; x hunter,n represents the position of the nth predator individual at the perceptible boundary; represents the variable of the rth predator in the dth dimension that the predator can perceive, r3 and r4 are two random numbers in the range [0,1];

[0055] S303: The predator adaptively tracks the prey and updates the position of the individual predator. The specific process is as follows:

[0056] The hunter group is randomly divided into two categories, ordinary hunters and elite hunters. There are differences in hunting ability between ordinary hunters and elite hunters. For the a-th ordinary hunter, a=p+1,…,N, where p∈[1,N], its position update process is shown as follows:

[0057]

[0058] in, represents the position of the ath common predator, Δda is the distance between the prey and the ath common predator; F is the symbol indicating the search direction, and its value is ±1; r5 is a random number in the range of [0,1];

[0059] For the b-th elite hunter, b = 1, ..., p, its position update process, that is, the active response to the prey escape, is as follows:

[0060]

[0061] in, represents the position of the bth elite hunter, F is the sign indicating the search direction, and its value is ±1; β is a constant indicating the proficiency of the hunter competition, r6, r7 and r8 are random numbers in the range of [0,1], is the competition intensity of the b-th elite hunter, Δd b represents the distance between the prey and the bth elite hunter;

[0062] S304: The prey is adaptively escaped from tracking and the position of the prey individual is updated. The specific process is as follows:

[0063] First, use the following formula to update the position of the prey in space:

[0064] x' prey =x prey +FβDI prey x prey +Fα·r9|1-cos2πr 10 |·Δd b |cos2πr 11 | (19)

[0065] Among them, x' prey is the updated prey position, r9, r 10 and r 11 are random numbers in the range [0,1];

[0066] For prey, it needs to escape from the hunter, that is, move to the safest place in space, that is, the place with the lowest fitness. This process is achieved through the BFGS quasi-Newton method; the BFGS quasi-Newton method uses the first-order derivative of the objective function and the approximate Hessian matrix information to achieve optimization in space through iterative solution;

[0067] S305, after updating the prey and predator groups, calculate and compare the fitness of all individuals, and take the individual with the best fitness as the prey individual to complete the update of all individuals in the entire population; return to S302, and continue to perform the t+1th iteration update of the identification parameter process based on the adaptive hunting game strategy;

[0068] When the parameter identification iterative process based on the adaptive hunting game strategy reaches the iterative termination condition, the parameter identification update iterative process is stopped, and the prey individual is output as the parameter identification result; wherein the iterative termination condition is to satisfy at least one of the following iterative update termination conditions: condition 1 is that the number of parameter update iterations is equal to the specified maximum number of iterative updates; condition 2 is that the fitness of the prey individual is less than the specified threshold;

[0069] S306, completing parameter identification and obtaining the final polarization characteristic model of the proton exchange membrane electrolyzer, as follows:

[0070] The parameter identification result obtained when the iteration is terminated is brought into the polarization characteristic model established in step S1, and finally a proton exchange membrane electrolyzer polarization characteristic model that can be used for industrial applications is obtained.

[0071] Furthermore, in step S304, the BFGS quasi-Newton method is used to use the first-order derivative of the objective function and the approximate Hessian matrix information to achieve optimization in space through iterative solution. The specific process is as follows:

[0072] The mean square error (MSE) between the measured output and the model output is used as the objective function of the quasi-Newton method, and its expression is:

[0073]

[0074] The first-order derivative of the objective function is:

[0075]

[0076] Among them, x d is the independent variable in the objective function, that is, the variable of each dimension, x d =x d ; Represents y i x d The partial derivative of It is denoted as g;

[0077] The approximate Hessian matrix is ​​obtained by formula (22):

[0078]

[0079] Where B is the approximation of the Hessian matrix, k represents the number of iterations of the BFGS method, the superscript T represents the transpose of the matrix, and s k With v k Obtained by the following formula:

[0080]

[0081] Among them, x prey,k represents the position of the prey in the kth iteration, sk is the search step vector, g k is the first-order derivative vector in the kth iteration calculated by equation (21);

[0082] In order to achieve a sufficient reduction in the objective function, x in formula (23) prey,k Using Armijo line search, we can determine the maximum value of w that satisfies equation (25) by looping, that is, we can determine the step size of the search and obtain x prey,k+1 ; The Armijo line search loop is shown in equations (24) to (26):

[0083]

[0084] x prey,k+1 =x prey,k +ρ w d k (25)

[0085]

[0086] Wherein, w represents the number of line search cycles. In the cycle, w is continuously increased to obtain the maximum w value that satisfies equation (24), that is, the step length of the search can be determined;

[0087] Among them, ρ and σ are constants, d k is the search direction vector;

[0088] The BFGS method can be cycled once to get a new x prey,k When the value satisfies the termination condition of the BFGS method, the iteration of the BFGS method is completed, and a new solution in the space is obtained, and the update of the prey position is completed; wherein, the termination condition of the BFGS method is to meet one of the following conditions: condition 1 is that the size of the first-order derivative is less than the specified threshold, and condition 2 is that the number of cycles of the BFGS method reaches the set maximum value;

[0089] The new solution obtained by the BFGS method is compared with the current optimal solution. If the new solution is better than the optimal solution, the new solution is used to replace the optimal solution as the new prey individual to complete the update of the prey individual.

[0090] Beneficial technical effects brought by the present invention:

[0091] Based on the existing polarization characteristic model of proton exchange membrane electrolyzer, the present invention combines the operation mechanism of the electrolyzer to improve the polarization characteristic model and enhance the modeling accuracy. At the same time, a parameter identification method based on an adaptive hunting game strategy is proposed. Under the premise of ensuring accurate parameter optimization, the rapid convergence of the parameter identification process is achieved, and the number of iterations required for parameter identification under small-scale populations is reduced. BRIEF DESCRIPTION OF THE DRAWINGS

[0092] Figure 1 A flow chart of a parameter identification method for improving the polarization characteristic model of a proton exchange membrane electrolyzer in the present invention; DETAILED DESCRIPTION

[0093] The specific implementation of the present invention is further described below in conjunction with specific embodiments:

[0094] A parameter identification method for improving the polarization characteristic model of proton exchange membrane electrolyzer, such as Figure 1 As shown, the following steps are included:

[0095] S1. Establish a polarization characteristic model of a proton exchange membrane electrolyzer and determine the model parameters to be identified;

[0096] The established polarization characteristic model of proton exchange membrane electrolyzer is as follows:

[0097] The proton exchange membrane electrolyzer is composed of several units, and the working voltage of each unit is u cell Expressed as open circuit voltage u ocv , activation overvoltage η act 、Ohmic overvoltage η ohm Concentration difference overvoltage η conc The sum of the four is shown in the following formula:

[0098] u cell =u ocv +η act +η ohm +η conc (1)

[0099] Open circuit voltage u ocv The expression is:

[0100]

[0101] Where T is the thermodynamic temperature scale, the unit is K, measured by experiment; R is the gas constant, generally taken as 8.3145 J / (mol·K); F is the Faraday constant, generally taken as 96485 C / mol; k wo is the correction coefficient, to be identified; p H2 is the hydrogen pressure in the electrolyzer, p O2 is the oxygen pressure in the electrolytic cell, and its expression is:

[0102] p H2 =P ca -p H2O (3)

[0103] p O2 =P an -p H2O (4)

[0104] Among them, P ca and P an are the absolute pressures measured at the hydrogen separator and the oxygen separator, respectively, p H2O is the water pressure in the electrolytic cell, and its expression is:

[0105]

[0106] Activation overvoltage η act The expression is:

[0107]

[0108] Among them, α an is the anode charge transfer coefficient, i is the current flowing through each electrolytic unit, i 0,an is the anode current density, and its expression is:

[0109]

[0110] Among them, i 0,an,std is the anode current density under standard conditions, E exc The activation energy required for electron transport to the anode;

[0111] Ohmic overvoltage η ohm The expression is:

[0112]

[0113] Among them, σ mem is the conductivity of the proton exchange membrane, and its expression is:

[0114]

[0115] Among them, σ mem,std is the conductivity of the proton exchange membrane under standard conditions, E pro The activation energy required for proton transport in the membrane;

[0116] Since the rated operating current density in the experiment is not high enough to produce concentration difference overvoltage η conc , so the value is 0.

[0117] The six model parameters that need to be identified in the polarization characteristic model of the proton exchange membrane electrolyzer are: correction coefficient k wo , anode charge transfer coefficient α an , the anode current density i under standard conditions 0,an,std , the activation energy E required for anode electron transfer exc , membrane resistivity σ under standard conditions mem,std and the activation energy E required for proton transport in the membranepro .

[0118] S2. Obtain a measurable operating variable data set of the proton exchange membrane electrolyzer, including a measurable input and output variable data set;

[0119] Specifically, the experimental data of measurable parameters of the electrolyzer and the corresponding working status of the electrolyzer are obtained and recorded through actual operating condition testing of the proton exchange membrane electrolyzer. The measurable parameters include experimental data such as voltage, current, temperature and pressure when the electrolyzer is operating.

[0120] S3. Based on the adaptive hunting game strategy, the parameters of the polarization characteristic model of the proton exchange membrane electrolyzer are identified to obtain the final polarization characteristic model of the proton exchange membrane electrolyzer.

[0121] Step S3 includes the following sub-steps:

[0122] S301, determine the reasonable range of the model parameters to be identified, and randomly generate initial population individuals within the range, specifically:

[0123] The range of the model parameters to be identified is given by existing public data or based on experience, and its expression is:

[0124]

[0125] Among them, lb1, lb2, lb3, lb4, lb5 and lb6 represent the lower bounds of the parameters to be identified, and ub1, ub2, ub3, ub4, ub5 and ub6 represent the upper bounds of the parameters to be identified;

[0126] According to the reasonable range of the model parameters to be identified, the initial population individuals are randomly generated, and the population number is set to N, that is, the number of individuals in the population. Each individual is a candidate solution, containing variables of 6 dimensions, representing the 6 parameters that need to be identified; the variable expression of the dth dimension of each candidate solution is:

[0127] x d =lb d +(ub d -lb d )·rand,d=1,2,...,6 (11)

[0128] Among them, x d represents the variable of the dth dimension, rand represents a random number in the range of [0,1], which is randomly determined when generating the variables of each dimension of all individuals;

[0129] S302. Evaluate individual danger intensity and competition intensity, specifically:

[0130] The root mean square error RMSE between the measured output and the model output is used as the fitness function f fitness (x k ), calculate the fitness of each individual, the expression is:

[0131]

[0132]

[0133] Among them, x k represents the position of the kth individual, x k =[x 1 ,x 2 ,x 3 ,x 4 ,x 5 ,x 6 ], h is the number of sample points contained in the measured data sequence, y i is the measured voltage output of the i-th sample point, is the model voltage output corresponding to the i-th sample point; the model voltage output can be calculated by the current, temperature, pressure and other data measured in step S2;

[0134] The smaller the fitness of an individual, the better it is. The individual with the smallest fitness is used as prey, and the other individuals in the population are used as predators.

[0135] For prey, DI prey It is defined as the intensity of danger perceived by the prey individual, and its expression is:

[0136]

[0137] Where m represents the number of predators that the prey can perceive, which is determined by roulette; x prey and x hunter,m Represent the positions of the prey and the mth predator at the perceptible boundary respectively; r1 and r2 are two random numbers in the range [0,1], and D represents the number of dimensions of the candidate solution, that is, the number of parameters. is the variable of the d-th dimension of the prey, is the variable of the jth predator in the dth dimension that the prey can perceive, α represents the endurance performance of the prey, which decreases with the increase of the number of iterations, and the expression is:

[0138]

[0139] Among them, C is the endurance constant, t is the current number of iterations, and t max The maximum number of iterations set;

[0140] For predators, DI hunter It is defined as the competition intensity perceived by individual predators and is expressed as:

[0141]

[0142] Where n represents the number of surrounding predators that the predator can perceive, which is determined by roulette; x hunter,n represents the position of the nth predator individual at the perceptible boundary; represents the variable of the rth predator in the dth dimension that the predator can perceive, r3 and r4 are two random numbers in the range [0,1];

[0143] S303: The predator adaptively tracks the prey and updates the position of the individual predator. The specific process is as follows:

[0144] The hunter group is randomly divided into two categories, ordinary hunters and elite hunters. There are differences in hunting ability between ordinary hunters and elite hunters. For the a-th ordinary hunter, a=p+1,…,N, where p∈[1,N], its position update process is shown as follows:

[0145]

[0146] in, represents the position of the ath common predator, Δd a is the distance between the prey and the ath common predator; F is the symbol indicating the search direction, and its value is ±1; r5 is a random number in the range of [0,1];

[0147] For the b-th elite hunter, b = 1, ..., p, its position update process, that is, the active response to the prey escape, is as follows:

[0148]

[0149] in, represents the position of the bth elite hunter, F is the sign indicating the search direction, and its value is ±1; β is a constant indicating the proficiency of the hunter competition, r6, r7 and r8 are random numbers in the range of [0,1], is the competition intensity of the b-th elite hunter, Δd b represents the distance between the prey and the bth elite hunter;

[0150] S304: The prey is adaptively escaped from tracking and the position of the prey individual is updated. The specific process is as follows:

[0151] First, use the following formula to update the position of the prey in space:

[0152] x' prey =x prey +FβDI prey x prey +Fα·r9|1-cos2πr 10 |·Δd b |cos2πr 11 | (19)

[0153] Among them, x' prey is the updated prey position, r9, r 10 and r 11 are random numbers in the range [0,1];

[0154] For prey, it needs to escape from the hunter, that is, move to the safest place in space, that is, the place with the lowest fitness. This process is achieved through the BFGS quasi-Newton method; the BFGS quasi-Newton method uses the first-order derivative of the objective function and the approximate Hessian matrix information to achieve optimization in space through iterative solution. The specific process is as follows:

[0155] The mean square error (MSE) between the measured output and the model output is used as the objective function of the quasi-Newton method, and its expression is:

[0156]

[0157] The first-order derivative of the objective function is:

[0158]

[0159] Among them, x d is the independent variable in the objective function, that is, the variable of each dimension, x d =x d ; Represents y i x d The partial derivative of It is denoted as g;

[0160] The approximate Hessian matrix is ​​obtained by formula (22):

[0161]

[0162] Where B is the approximation of the Hessian matrix, k represents the number of iterations of the BFGS method, the superscript T represents the transpose of the matrix, and s k With v k Obtained by the following formula:

[0163]

[0164] Among them, x prey,k represents the position of the prey in the kth iteration, sk is the search step vector, g k is the first-order derivative vector in the kth iteration calculated by equation (21);

[0165] In order to achieve a sufficient reduction in the objective function, x in formula (23) prey,k Using Armijo line search, we can determine the maximum value of w that satisfies equation (25) by looping, that is, we can determine the step size of the search and obtain x prey,k+1 ; The Armijo line search loop is shown in equations (24) to (26):

[0166]

[0167] x prey,k+1 =x prey,k +ρ w d k (25)

[0168]

[0169] Wherein, w represents the number of line search cycles. In the cycle, w is continuously increased to obtain the maximum w value that satisfies equation (24), that is, the step length of the search can be determined;

[0170] Among them, ρ and σ are constants, which are generally determined after testing based on specific problems. k is the search direction vector;

[0171] The BFGS method can be cycled once to get a new x prey,k When the value satisfies the termination condition of the BFGS method, the iteration of the BFGS method is completed, and a new solution in the space is obtained, and the update of the prey position is completed; wherein, the termination condition of the BFGS method is to meet one of the following conditions: condition 1 is that the size of the first-order derivative is less than the specified threshold, and condition 2 is that the number of cycles of the BFGS method reaches the set maximum value;

[0172] The new solution obtained by the BFGS method is compared with the current optimal solution. If the new solution is better than the optimal solution, the new solution is used to replace the optimal solution as the new prey individual to complete the update of the prey individual.

[0173] S305, after updating the prey and predator groups, calculate and compare the fitness of all individuals, and take the individual with the best fitness as the prey individual to complete the update of all individuals in the entire population; return to S302, and continue to perform the t+1th iteration update of the identification parameter process based on the adaptive hunting game strategy;

[0174] When the parameter identification iterative process based on the adaptive hunting game strategy reaches the iterative termination condition, the parameter identification update iterative process is stopped, and the prey individual is output as the parameter identification result; wherein the iterative termination condition is to satisfy at least one of the following iterative update termination conditions: condition 1 is that the number of parameter update iterations is equal to the specified maximum number of iterative updates; condition 2 is that the fitness of the prey individual is less than the specified threshold;

[0175] S306, completing parameter identification and obtaining the final polarization characteristic model of the proton exchange membrane electrolyzer, as follows:

[0176] The parameter identification result obtained when the iteration is terminated is brought into the polarization characteristic model established in step S1, and finally a proton exchange membrane electrolyzer polarization characteristic model that can be used for industrial applications is obtained.

[0177] Of course, the above description is not a limitation of the present invention, and the present invention is not limited to the above examples. Changes, modifications, additions or substitutions made by technicians in this technical field within the essential scope of the present invention should also fall within the protection scope of the present invention.

Claims

1. A parameter identification method for improving the polarization characteristic model of a proton exchange membrane electrolyzer, characterized in that: The following steps are involved: S1. Establish a polarization characteristic model of a proton exchange membrane electrolyzer and determine the model parameters to be identified; S2. Obtain a measurable operating variable data set of the proton exchange membrane electrolyzer, including a measurable input and output variable data set; S3. Based on the adaptive hunting game strategy, the parameters of the polarization characteristic model of the proton exchange membrane electrolyzer are identified, and then the final polarization characteristic model of the proton exchange membrane electrolyzer is obtained; In step S1, the established proton exchange membrane electrolyzer polarization characteristic model is specifically: The proton exchange membrane electrolyzer is composed of several units, and the working voltage of each unit is u cell Expressed as open circuit voltage u ocv , activation overvoltage η act 、Ohmic overvoltage η ohm Concentration difference overvoltage η conc The sum of the four is shown in the following formula: you cell =u ocv +n act +n ohm +n conc (1) Open circuit voltage u ocv The expression is: Where T is the thermodynamic temperature scale, the unit is K, measured by experiment; R is the gas constant; F is the Faraday constant; k wo is the correction factor; p H2 is the hydrogen pressure in the electrolyzer, p O2 is the oxygen pressure in the electrolytic cell, and its expression is: p H2 =P ca -p H2O (3) p O2 =P an -p H2O (4) Among them, P ca and P an are the absolute pressures measured at the hydrogen separator and the oxygen separator, respectively, p H2O is the water pressure in the electrolytic cell, and its expression is: Activation overvoltage η act The expression is: Among them, α an is the anode charge transfer coefficient, i is the current flowing through each electrolytic unit, i 0,an is the anode current density, and its expression is: Among them, i 0,an,std is the anode current density under standard conditions, E exc The activation energy required for electron transport to the anode; Ohmic overvoltage η ohm The expression is: Among them, σ mem is the conductivity of the proton exchange membrane, and its expression is: Among them, σ mem,std is the conductivity of the proton exchange membrane under standard conditions, E pro The activation energy required for proton transport in the membrane; Since the rated operating current density in the experiment is not high enough to produce concentration difference overvoltage η conc , so the value is 0; The six model parameters that need to be identified in the polarization characteristic model of the proton exchange membrane electrolyzer are: correction coefficient k wo , anode charge transfer coefficient α an , the anode current density i under standard conditions 0,an,std , the activation energy E required for anode electron transfer exc , membrane resistivity σ under standard conditions mem,std and the activation energy E required for proton transport in the membrane pro ; The step S3 comprises the following sub-steps: S301, determine the reasonable range of the model parameters to be identified, and randomly generate initial population individuals within the range, specifically: The range of the model parameters to be identified is given by existing public data or based on experience, and its expression is: Among them, lb1, lb2, lb3, lb4, lb5 and lb6 represent the lower bounds of the parameters to be identified, and ub1, ub2, ub3, ub4, ub5 and ub6 represent the upper bounds of the parameters to be identified; According to the reasonable range of the model parameters to be identified, the initial population individuals are randomly generated, and the population number is set to N, that is, the number of individuals in the population. Each individual is a candidate solution, containing variables of 6 dimensions, representing the 6 parameters that need to be identified; the variable expression of the dth dimension of each candidate solution is: x d =lb d +(ub d -lb d )·rand,d=1,2,...,6 (11) Among them, x d represents the variable of the dth dimension, rand represents a random number in the range of [0,1], which is randomly determined when generating the variables of each dimension of all individuals; S302. Evaluate individual danger intensity and competition intensity, specifically: The root mean square error RMSE between the measured output and the model output is used as the fitness function f fitness (x k ), calculate the fitness of each individual, the expression is: Among them, x k represents the position of the kth individual, x k =[x 1 ,x 2 ,x 3 ,x 4 ,x 5 ,x 6 ], h is the number of sample points contained in the measured data sequence, y i is the measured voltage output of the i-th sample point, is the model voltage output corresponding to the i-th sample point; The individual with the lowest fitness is regarded as prey, and the other individuals in the population are regarded as predators; For prey, DI prey It is defined as the intensity of danger perceived by the prey individual, and its expression is: Where m represents the number of predators that the prey can perceive, which is determined by roulette; x prey and x hunter,m Represent the positions of the prey and the mth predator at the perceptible boundary respectively; r1 and r2 are two random numbers in the range [0,1], and D represents the number of dimensions of the candidate solution, that is, the number of parameters. is the variable of the d-th dimension of the prey, is the variable of the jth predator in the dth dimension that the prey can perceive, α represents the endurance performance of the prey, which decreases with the increase of the number of iterations, and the expression is: Among them, C is the endurance constant, t is the current number of iterations, and t max The maximum number of iterations set; For predators, DI hunter It is defined as the competition intensity perceived by individual predators and is expressed as: Where n represents the number of surrounding predators that the predator can perceive, which is determined by roulette; x hunter,n represents the position of the nth predator individual at the perceptible boundary; represents the variable of the rth predator in the dth dimension that the predator can perceive, r3 and r4 are two random numbers in the range [0,1]; S303: The predator adaptively tracks the prey and updates the position of the individual predator. The specific process is as follows: The hunter group is randomly divided into two categories, ordinary hunters and elite hunters. There are differences in hunting ability between ordinary hunters and elite hunters. For the a-th ordinary hunter, a=p+1,…,N, where p∈[1,N], its position update process is shown as follows: in, represents the position of the ath common predator, Δd a is the distance between the prey and the ath common predator; F is the symbol indicating the search direction, and its value is ±1; r5 is a random number in the range of [0,1]; For the b-th elite hunter, b = 1, ..., p, its position update process, that is, the active response to the prey escape, is as follows: in, represents the position of the bth elite hunter, F is the sign indicating the search direction, and its value is ±1; β is a constant indicating the proficiency of the hunter competition, r6, r7 and r8 are random numbers in the range of [0,1], is the competition intensity of the b-th elite hunter, Δd b represents the distance between the prey and the bth elite hunter; S304: Prey adaptively escapes tracking and updates the location of the prey individual. The specific process is as follows: First, use the following formula to update the position of the prey in space: x' prey =x prey +FβDI prey x prey +Fα·r9|1-cos2πr 10 |·Δd b |cos2πr 11 | (19) Among them, x' prey is the updated prey position, r9, r 10 and r 11 are random numbers in the range [0,1]; For prey, it needs to escape from the hunter, that is, move to the safest place in space, that is, the place with the lowest fitness. This process is achieved through the BFGS quasi-Newton method; the BFGS quasi-Newton method uses the first-order derivative of the objective function and the approximate Hessian matrix information to achieve optimization in space through iterative solution; S305, after updating the prey and predator groups, calculate and compare the fitness of all individuals, and take the individual with the best fitness as the prey individual to complete the update of all individuals in the entire population; return to S302, and continue to perform the t+1th iteration update of the identification parameter process based on the adaptive hunting game strategy; When the parameter identification iterative process based on the adaptive hunting game strategy reaches the iterative termination condition, the parameter identification update iterative process is stopped, and the prey individual is output as the parameter identification result; wherein the iterative termination condition is to satisfy at least one of the following iterative update termination conditions: condition 1 is that the number of parameter update iterations is equal to the specified maximum number of iterative updates; condition 2 is that the fitness of the prey individual is less than the specified threshold; S306, completing parameter identification and obtaining the final polarization characteristic model of the proton exchange membrane electrolyzer, as follows: The parameter identification result obtained when the iteration ends is brought into the polarization characteristic model established in step S1, thereby finally obtaining a polarization characteristic model of a proton exchange membrane electrolyzer that can be used for industrial applications; In step S304, the BFGS quasi-Newton method is used to use the first-order derivative of the objective function and the approximate Hessian matrix information to achieve optimization in space through iterative solution. The specific process is as follows: The mean square error (MSE) between the measured output and the model output is used as the objective function of the quasi-Newton method, and its expression is: The first-order derivative of the objective function is: Among them, x d is the independent variable in the objective function, that is, the variable of each dimension, x d =x d ; Represents y i x d The partial derivative of It is denoted as g; The approximate Hessian matrix is ​​obtained by formula (22): Where B is the approximation of the Hessian matrix, k represents the number of iterations of the BFGS method, T represents the transpose of the matrix, and s k With v k Obtained by the following formula: Among them, x prey,k represents the position of the prey in the kth iteration, s k is the search step vector, g k is the first-order derivative vector in the kth iteration calculated by equation (21); In order to achieve a sufficient reduction in the objective function, x in formula (23) prey,k Using Armijo line search, we can determine the maximum value of w that satisfies equation (25) by looping, that is, we can determine the step size of the search and obtain x prey,k+1 ; The Armijo line search loop is shown in equations (24) to (26): x prey,k+1 =x prey,k +r w d k (25) Wherein, w represents the number of line search cycles. In the cycle, w is continuously increased to obtain the maximum w value that satisfies equation (24), that is, the step length of the search can be determined; Among them, ρ and σ are constants, d k is the search direction vector; The BFGS method can be cycled once to get a new x prey,k When the value satisfies the termination condition of the BFGS method, the iteration of the BFGS method is completed, and a new solution in the space is obtained, and the update of the prey position is completed; wherein, the termination condition of the BFGS method is to meet one of the following conditions: condition 1 is that the size of the first-order derivative is less than the specified threshold, and condition 2 is that the number of cycles of the BFGS method reaches the set maximum value; The new solution obtained by the BFGS method is compared with the current optimal solution. If the new solution is better than the optimal solution, the new solution is used to replace the optimal solution as the new prey individual to complete the update of the prey individual.

2. The method for parameter identification of an improved proton exchange membrane electrolyzer polarization characteristic model according to claim 1, characterized in that: In step S2, experimental data of measurable parameters of the electrolyzer and corresponding working status of the electrolyzer are obtained and recorded through actual working condition testing of the proton exchange membrane electrolyzer, and the measurable parameters include voltage, current, temperature and pressure when the electrolyzer is in operation.

Images