Three-level power generation system multi-step prediction control method based on rolling time domain estimation

Abstract

The invention discloses a three-level power generation system multi-step prediction control method based on rolling time domain estimation. The method comprises the following steps: establishing a mathematical model under a new topology based on a three-level structure adopted by a machine-side converter and a network-side converter; discretizing the model; sampling current, voltage and speed and performing coordinate transformation; designing an outer ring model prediction speed controller; designing an inner ring model prediction current controller; selecting an optimal voltage vector enabling the cost function value to be minimum according to the constructed cost function; the dq-axis observer is used for estimating dq-axis disturbance information of the system, so that the designed current loop controller has relatively strong suppression on parameter mismatching; estimating state information of the system and disturbance information involved in each step of multi-step prediction by means of a differential GPIO observer, and constructing a new state space model; according to the method, the problem of large calculation amount is solved, and the anti-interference performance of the system is effectively enhanced.

Description

Technical Field

[0001] The present invention relates to the technical field of wind power generation, and in particular to a multi-step predictive control method for a three-level power generation system based on rolling time domain estimation. Background Art

[0002] With the rapid development of industrialization, environmental pollution and energy shortages have become increasingly prominent. As an efficient, clean, pollution-free, and single-management energy form, how to use natural energy to generate electricity efficiently has become a research hotspot that has attracted widespread attention from all countries.

[0003] In recent years, permanent magnet direct drive power generation system has received extensive attention in the field of wind power generation. This system has the advantages of high power generation efficiency at low wind speed and direct drive (without gearbox). However, the traditional PID control method requires accurate modeling of the wind power generation system, and the complex characteristics of the wind power generation system make accurate modeling difficult, resulting in poor control effect. In response to these problems, domestic and foreign scholars have proposed a variety of advanced control strategies, such as predictive current control, sliding mode variable structure control, and anti-disturbance control.

[0004] Compared with PID control, the controller design of finite set model predictive control is more flexible. The cost function can be constructed according to the actual control target. On the other hand, finite set model predictive control can directly act on the optimal switching state of the converter according to the constructed cost function output without modulation, which greatly speeds up the dynamic response of the system. Because the self-disturbance rejection control adopts a two-degree-of-freedom structure, it can well balance the tracking performance and anti-disturbance performance. The dq-axis differential observer is used to estimate the external disturbance and uncertainty factors of the system in real time, and compensate them on the controller side, which enhances the disturbance suppression ability of the system while meeting the tracking control performance.

[0005] In order to improve the dynamic response and anti-disturbance performance of the permanent magnet direct-drive wind power generation system, under the three-level topology, predictive control is used in the outer loop for prediction and estimation; the inner loop adopts multi-step predictive control based on spherical decoding, and the corresponding dq axis differential observer is designed for compensation. The traditional anti-disturbance control can only estimate and compensate the constant mismatch disturbance well, but cannot estimate and compensate the time-varying disturbance of the mismatch disturbance well. Therefore, the dq axis differential observer is designed to suppress the disturbance, so that the designed current loop controller has a strong suppression of parameter mismatch. On the other hand, the traditional two-level topology can only reduce the output harmonics by increasing the switching frequency, but too high a switching frequency will increase the switching loss, while the multi-level topology can make the switching device switch once in each cycle to achieve the same effect as the traditional converter switch several times, so at the same switching frequency, more harmonic components can be eliminated. Therefore, in order to reduce the harmonic components, the machine-side and grid-side converter of the permanent magnet direct-drive wind power generation system adopts a diode clamped three-level topology, which improves the problems of large switching loss and high harmonics to a certain extent. Summary of the invention

[0006] In order to overcome the shortcomings of the prior art, the present invention provides a multi-step predictive control method for a three-level power generation system based on rolling time domain estimation. Under the three-level topology, the outer loop adopts predictive control for predictive estimation; the inner loop adopts multi-step predictive control based on spherical decoding, and designs a corresponding dq axis differential observer for compensation, thereby solving the technical problem.

[0007] In order to solve the above technical problems, the technical solution provided by the present invention is as follows:

[0008] A multi-step predictive control method for a three-level power generation system based on a rolling horizon estimation comprises the following steps:

[0009] Step 1: Determine the given speed value ω of the machine side speed loop by capturing the wind speed ref ;

[0010] Step 2: Establish the machine side mathematical model;

[0011] Step 3: Sample the generator side current and speed and perform Clark-Park transformation. The process is as follows:

[0012] First, the three-phase current i of the machine-side converter is collected in real time a ,i b ,i c , and then the collected i a ,i b ,i c Perform Clark transformation to obtain the component value i in the α-β coordinate system α ,i β, and perform Park transformation in the α-β coordinate system, the component value i in the α-β coordinate system can be α ,i β Convert to component value i in dq coordinate system d ,i q ;

[0013] Step 4: Establish a discrete permanent magnet synchronous generator speed prediction model. The process is as follows:

[0014] The dynamic mathematical equation of the permanent magnet synchronous generator is:

[0015]

[0016] in P n is the number of pole pairs of the surface-mounted permanent magnet synchronous generator, ψ f is the flux linkage of the surface mounted permanent magnet synchronous generator, T m is the torque of the surface mounted permanent magnet synchronous generator, J is the moment of inertia, and B is the motor friction coefficient.

[0017] right Perform Laplace transform to get the first-order model, and then perform Euler discretization to get the discrete-time model:

[0018] ω r (k) = ai q (k-1)-bω r (k-1)

[0019] ω r (k+1)=ai q (k)-bω r (k)

[0020] where ω r (k) is the actual speed of the motor at time kT, b=-e -TB / J , the speed prediction model is obtained by rearranging the above formula as follows:

[0021]

[0022] Among them, ω m (k+1) is the predicted output at time (k+1)T, is the q-axis current control increment of the schedule permanent magnet synchronous generator at time kT.

[0023] When the prediction step is extended to N steps, we can get:

[0024]

[0025] Among them, A=-b, B=a; Δiq (k) is the q-axis current control increment of the schedule-mounted permanent magnet synchronous generator at time kT, A i Represents the matrix obtained by raising the matrix A to the i-th power.

[0026] can be simplified to:

[0027] W ωm (k) = W ωr (k)+GΔI q (k)

[0028] in,

[0029] Step 5: Establish a discrete permanent magnet synchronous generator current prediction model. The process is as follows:

[0030] According to the coordinate transformation theory, the voltage equation of the surface-mount PMSG mathematical model in the dq coordinate system can be obtained as follows:

[0031]

[0032] Where: u d ,u q The d and q axis voltage components of the stator voltage corresponding to 27 switching states; i d ,i q represents the d and q axis components of the stator current in the dq coordinate system; L s is the stator inductance of the surface-mounted permanent magnet synchronous generator in the dq coordinate system; R s represents the stator resistance; ψ f represents the permanent magnet flux; ω re Represents electrical angular velocity.

[0033] The forward Euler formula is used to discretize the current state equation and after one step of delay, the discrete permanent magnet synchronous generator current prediction model is obtained as follows:

[0034]

[0035] Among them, i d (k+1), i q (k+1) is the d,q axis component of the stator current in the dq coordinate system after delay compensation; i d (k+2), i q (k+2) are the predicted values ​​of the d and q axis stator currents at time k+1 respectively; ω re (k+1) is approximately equal to ω re (k);

[0036] Taking into account the uncertainty of system parameters and the influence of external disturbances, the above current equation discretized by Euler is transformed into a vector form:

[0037] i dq (k+1)=Ai dq (k)+BPu(k)+C+B d d(k)

[0038] Among them, among them,

[0039] i dq =[i d i q ] T ,u(k)=[S a S b S c ],

[0040] d(k)=[d d (k) d q (k)] T ,

[0041] Where S a ,S b ,S c Respectively represent the states of the upper and lower switches of the three bridge arms of the inverter, u dc Represents the voltage across the DC bus capacitor.

[0042] When the prediction step is extended to N steps, it is assumed that the rotor speed of the motor remains unchanged within the prediction range, and it is defined as:

[0043]

[0044] Where I(k) represents the matrix of the predicted dq-axis current values ​​for the next N time steps at time KT, I * (k) represents the dq-axis reference current matrix for the next N time steps at time KT, A matrix representing the lumped perturbation estimate at time KT for N time steps into the future.

[0045] Thus, we can get:

[0046]

[0047] in

[0048] Step 6, construct the cost function;

[0049] A cost function is constructed to better represent the tracking effect of the predicted value on the expected value. The cost function J1 is as follows:

[0050]

[0051] in, Indicates the reference value of the stator current d, q axis components; i dq (l+1) are the predicted values ​​of d-axis and q-axis stator currents at time l+1, u is the predicted switching frequency, and λ represents the weight coefficient.

[0052] In order to comply with the calculation rules of the sphere decoding algorithm, it is transformed into the following J1 form:

[0053]

[0054] in, Bundle Substitute the cost function J1 and further sort J1 using the matching method to obtain:

[0055] J1=U(k) T HU(k)+2Θ(k) T U(k)+θ(k)=(U(k)+H -1 Θ(k) T H(U(k)+H -1 Θ(k))+const(k) where, const(k) is a constant that is independent of U(k), so the constant term is independent of the switch switching sequence. At this time, the optimization problem can be transformed into:

[0056]

[0057] Since the matrix H is a positive definite symmetric matrix, the equation can be solved to obtain U unc (k)=-H -1 Θ(k), we get the new cost function J'1:

[0058] J'1=(U(k)-U unc (k)) T H(U(k)-U unc (k))

[0059] Among them, U unc (k) is the optimal solution without constraints.

[0060] With the matrix H being a positive definite symmetric matrix, it can be concluded that its inverse matrix is ​​also a positive definite matrix; therefore, there exists a unique reversible lower triangular matrix V that satisfies:

[0061] H-1 =VV T

[0062] make Then the cost function J'1 is transformed into the following formula:

[0063]

[0064] Among them, V -1 is a lower triangular matrix.

[0065] A cost function is constructed to better represent the tracking effect of the predicted value on the expected value. The cost function J2 is as follows:

[0066]

[0067] It is necessary to predict the error in the future optimization time domain and compensate for it.

[0068] e(k)=ω r (k)-ω m (k) and e(k)=e(k+1)=…=e(k+n)

[0069] In order to smoothly transition to the set point, a first-order exponential change is adopted, and the reference trajectory equation is as follows:

[0070] y p (k+j|k)=ω ref (k+j)-α j [ω ref (k)-ω r (k)],j=1…N

[0071] Among them, ω ref (k) is the setting value, α is the weight coefficient, y p (k+1) is the expected output determined by the reference trajectory.

[0072] After prediction correction

[0073] ω' r (k+j|k)=y p (k+j|k)+e(k), j=1…N

[0074] Among them, ω' r (k+j|k) is the corrected reference speed, and e(k) represents the prediction error.

[0075] The above cost function is simplified to:

[0076]

[0077] in,

[0078] Step 7, selecting the optimal voltage vector;

[0079] First, the output voltage vector of the three-level inverter is determined by the switching state of the three bridge arms; then, the prediction value at the current moment can be obtained under the action of the prediction model; finally, the optimal voltage vector u is selected according to the designed cost function. opt_1 :

[0080] u opt_1 =argmin J'1

[0081] First, the reference speed at the current moment is collected; then, under the action of the prediction model, the speed prediction value at the current moment can be obtained; finally, the optimal voltage vector u is selected according to the designed cost function. opt_2 :

[0082] u opt_2 =argmin J2

[0083] Step 8, the dq axis differential observer process is as follows:

[0084] Taking into account the uncertainty of system parameters and the influence of external disturbances, the permanent magnet synchronous generator equation in step 5 can be sorted out as follows:

[0085] Let x1(k)=i d (k), then

[0086]

[0087] Among them, d d (k) represents the lumped disturbance, f1(k), f2(k), and f3(k) represent d d (k) first derivative, second derivative, third derivative.

[0088] make Design the corresponding observer, then

[0089]

[0090] Select appropriate observer gains β1, β2, β3, and β4 to achieve asymptotic convergence of the designed dq-axis differential observer. At this time, it is necessary to extract the disturbance information involved in each step of the multi-step FCS-MPC from the observer's estimated value. The specific extraction formula is as follows:

[0091]

[0092] Similarly, a similar observer is designed for the q-axis. Based on the disturbance information extracted from the dq-axis difference observer, the required spherical decoding prediction model can be easily obtained.

[0093] Step 9: Introduce state variable d ωl , determine the new state space model;

[0094] In the actual speed regulation process of the motor, the steady-state performance of the speed outer loop will be affected by external disturbances and uncertainty of model parameters. The mechanical motion equation is:

[0095]

[0096] Due to the uncertainty of model parameters and the influence of external disturbances, the variable d is introduced ωl , transforming the above mechanical motion equation into this form:

[0097]

[0098] in, represents the lumped disturbance on the generator side; It represents the control quantity of the machine-side controller after disturbance compensation; ω represents the mechanical angular velocity.

[0099] Let x 11 =ω,x 12 =d ωl is the new state variable, then the new state space model is expressed as:

[0100]

[0101] Among them, x 11 =ω represents the mechanical angular velocity, x 12 =d ωl represents the lumped disturbance, x 13 For x 12 The first derivative of , h1 represents d ωl The second-order differential of It represents the control quantity of the machine-side controller after disturbance compensation.

[0102] Step 10, differential GPIO observer design, the process is as follows:

[0103] According to the new state space model in step 9, the differential GPIO observer is designed in the form of:

[0104]

[0105] in, represents the error variable, represents the estimated value of ω, represents the estimated value of the lumped disturbance, is the estimated value of the first-order derivative of the lumped disturbance; It represents the control quantity of the machine-side controller after disturbance compensation;

[0106] β 11 ,β 12 ,β 13 ,β 14 represents the gain of the extended state observer;

[0107] make The forward Euler formula is used to discretize the above formula to obtain:

[0108]

[0109] At this point, we need to extract β from the estimated value of the differential GPIO observer 11 ,β 12 ,β 13 ,β 14 The perturbation information involved in each step of multi-step speed prediction. The specific extraction formula is as follows:

[0110]

[0111] Step 11, establishing a DC link mathematical model;

[0112] Step 12, establishing a network-side mathematical model;

[0113] Step 13, sampling and coordinate transformation of grid-side current and voltage;

[0114] Step 14, establishing a discrete inner loop power prediction model;

[0115] Step 15, constructing a cost function;

[0116] Step 16, selecting the optimal voltage vector;

[0117] Step 17, introduce state variable d ul , determine the new state space model;

[0118] Step 18, designing an extended state observer;

[0119] Step 19: Design of the grid-side outer loop control law.

[0120] Furthermore, in step 11, the current at the DC side capacitor nodes P, O, and N is expressed as:

[0121]

[0122] i c1 =i pm -i pg

[0123] i c1 +iom =i c2 +i og

[0124] i c2 +i nm =i ng

[0125] Where, C1, C2 represent DC filter capacitors; u c1 ,u c2 Represents the voltage on the DC bus capacitor; i c1 ,i c2 Represents the current flowing through the DC filter capacitor; i pm ,i om ,i nm Indicates the current flowing through nodes P, O, and N on the machine side, i pg ,i og ,i ng Indicates the current flowing to the grid-side nodes P, O, and N.

[0126] Furthermore, in step 12, the network-side mathematical model in the dq coordinate system is:

[0127]

[0128] Among them, u d ,u q is the component of the output voltage of the three-level inverter in the dq coordinate system; e d ,e q is the grid-side voltage component in the dq coordinate system; i d ,i q is the grid-side current component in the dq coordinate system; L represents the grid-side filter inductance; R represents the output equivalent resistance; ω ge Represents the grid angular velocity.

[0129] Furthermore, in step 14, the grid-side inverter adopts a voltage-oriented control method, so the grid-side inverter current equation based on grid voltage vector orientation can be expressed as:

[0130]

[0131] Among them, e d is the d-axis component of the grid-side voltage.

[0132] Furthermore, in step 15, a cost function is constructed, and the cost function J3 is as follows:

[0133]

[0134] Among them, P * ,Q *Indicates the reference value of active power and reactive power; P(k+1), Q(k+1) are (k+1)T s Predicted values ​​of active power and reactive power at each moment.

[0135] Furthermore, in step 16, a voltage vector u that minimizes the cost function J2 is selected from the 27 voltage vectors output by the grid-side inverter. opt_3 :

[0136] u opt_3 =argmin J3

[0137] Furthermore, in step 17, the state variable d is introduced ul , construct a new state space model;

[0138] Without considering the converter loss, the output power P of the machine-side rectifier is m It can be expressed as:

[0139] P m =u dc i m

[0140] Among them, u dc Represents the DC bus voltage, which can be expressed as u dc =u c1 +u c2 ;i m Indicates the current output from the machine-side converter to the DC bus;

[0141] The current flowing through the DC capacitor is:

[0142]

[0143] Wherein, C represents the DC side capacitance, which can be expressed as C=C1=C2; i g Indicates the current input to the grid-side inverter;

[0144] The active power P input from the DC side to the grid-side inverter is:

[0145] P=u dc i g

[0146] From the above equation we can get:

[0147]

[0148] Equivalent to

[0149]

[0150] in, represents the total disturbance of the grid-side voltage loop; b u0 It's about b u The estimated value of P * Indicates the active power reference value;

[0151] Let z1=u dc ; z2=d ul , then the new state space model is:

[0152]

[0153] Among them, h2 represents d ul The differential of u0 It's about b u The estimated value of P * Indicates the active power reference value.

[0154] Furthermore, in step 18, the improved extended state observer of the outer loop is expressed as follows:

[0155]

[0156] in, Indicates u dc An estimated value of represents the lumped disturbance d ul The estimated value of; l1, l2, l3 represent the gains of the improved extended state observer;

[0157] Furthermore, in step 19, the outer loop control law is designed as follows:

[0158] By selecting appropriate observer gains l1, l2, l3, the estimated value of the DC bus voltage can be obtained by the extended state observer designed in step 18. and an estimate of the outer loop lumped disturbance The estimated value obtained by the extended state state observer can be used for controller design, and the specific form is as follows:

[0159]

[0160] in, Indicates the reference value of the voltage outer loop; u u0 Indicates the grid-side controller output; k up represents the controller gain.

[0161] Design concept of the present invention:

[0162] The present invention adopts a three-level structure for the machine-side and grid-side converters, establishes a mathematical model under a new topology, discretizes the model, samples the current, voltage and speed, and performs coordinate transformation, designs an outer-loop model prediction speed controller, designs an inner-loop model prediction current controller, selects an optimal voltage vector that minimizes the cost function value according to the constructed cost function, estimates the information of the system dq-axis disturbance with the help of a dq-axis observer, so that the designed current loop controller has a strong suppression of parameter mismatch, estimates the system state information and the disturbance information involved in each step of multi-step prediction with the help of a differential GPIO observer, and constructs a new state space model. The present invention provides a spherical decoding multi-step prediction control method for a three-level power generation system based on rolling optimization, optimizes the traditional finite set model prediction control strategy, uses prediction control for the outer loop for prediction estimation, greatly reduces the problem of large calculation amount and effectively enhances the anti-disturbance performance of the system.

[0163] The beneficial effects of the present invention are:

[0164] The present invention adopts predictive control for the outer loop of the three-level power generation system to perform predictive estimation; the current loop adopts multi-step spherical decoding model current predictive control to replace the traditional PI controller, which has the advantage of not requiring a modulator, thereby improving the dynamic response speed of the three-level power generation system and greatly reducing the amount of calculation compared to the traditional FCS-MPC. At the same time, a dq axis differential observer is designed so that the designed current loop controller has a strong suppression of parameter mismatch. BRIEF DESCRIPTION OF THE DRAWINGS

[0165] Figure 1 It is the overall structure diagram of the three-level permanent magnet direct-drive wind power generation system;

[0166] Figure 2 is the space vector diagram of the three-level inverter;

[0167] Figure 3 This is the block diagram of the finite set model predictive current control based on the first-order active disturbance rejection control on the generator side;

[0168] Figure 4 It is a block diagram of the finite set model predictive power control based on the first-order active disturbance rejection control on the grid side;

[0169] Figure 5 This is the simulation diagram of the speed waveform when the wind speed increases at 0.6s;

[0170] Figure 6 This is the simulation diagram of the electromagnetic torque waveform when the wind speed increases at 0.6s;

[0171] Figure 7 This is the simulation diagram of the tracking waveform of the q-axis stator current component of the current loop on the mechanical side when the wind speed increases at 0.6s;

[0172] Figure 8 This is the voltage waveform simulation diagram of the voltage loop when the wind speed increases at 0.6s;

[0173] Fig. 9 This is the waveform simulation diagram of grid-side active power tracking when wind speed increases at 0.6s;

[0174] Fig.10 This is the voltage loop voltage waveform simulation diagram when the grid voltage changes at 0.8s;

[0175] Fig.11 This is a simulation diagram of the speed waveform when a disturbance is applied for 0.9 seconds under constant wind speed conditions. DETAILED DESCRIPTION

[0176] In order to make the technical solution of the present invention clearer, a detailed description is given below in conjunction with the accompanying drawings. The specific implementation examples described here are only used to explain the present invention and are not used to limit the present invention.

[0177] Reference Figures 1 to 11 , a multi-step predictive control method for a three-level power generation system based on a rolling horizon estimation, comprising the following steps:

[0178] Step 1: Determine the given speed value ω of the machine side speed loop ref ;

[0179] Step 2: Establish the machine-side mathematical model. The process is as follows:

[0180] The mathematical model of the permanent magnet synchronous generator in the dq coordinate system can be expressed as:

[0181]

[0182] Where: u d ,u q represents the dq-axis component of the stator voltage; i d ,i q Represents the dq-axis component of the stator current; L s is the stator inductance in the dq coordinate system of the surface-mounted permanent magnet synchronous generator, satisfying L s =L d =L q ; R s represents the stator resistance; ω re represents the electrical angular velocity; ψ f represents the permanent magnet flux.

[0183] The electromagnetic torque equation is:

[0184]

[0185] Among them, p nrepresents the number of pole pairs; T e Represents electromagnetic torque.

[0186] The mechanical motion equation is:

[0187]

[0188] Where, ω represents the mechanical angular velocity; J represents the moment of inertia; B represents the friction coefficient; T m Indicates the driving torque.

[0189] Step 3, sampling the machine side current and speed and performing coordinate transformation;

[0190] In order to achieve effective control of the permanent magnet synchronous generator on the machine side, double closed-loop control is adopted. The control of the current loop requires knowing the current information of the permanent magnet synchronous generator. Since the current information collected in real time is in the three-phase stationary coordinate system, in order to facilitate control, the current information in the three-phase stationary coordinate system needs to be converted to the dq coordinate system.

[0191] Clark Transformation:

[0192]

[0193] Park Transform:

[0194]

[0195] Among them, θ re is the electrical angle of rotation, satisfying

[0196] Step 4: Establish a discrete permanent magnet synchronous generator speed prediction model. The process is as follows:

[0197] 4.1: Determination of output voltage vector of three-level inverter:

[0198] Assume the three-phase sinusoidal voltage expression is:

[0199]

[0200] The inverter output voltage is defined as:

[0201]

[0202] but

[0203]

[0204] And U aN +U bN +U cN =0, then

[0205]

[0206] Among them, S a , S b , S c Corresponding to the three switching states of the three bridge arms of the inverter, U aN , U bN , U cN They correspond to the voltage between the inverter output phase A, phase B, phase C and the neutral point N, U oN Represents the inverter output voltage, U dc Represents the voltage across the DC bus capacitor.

[0207] The relationship between the switching states of the three bridge arms of the three-level inverter and the inverter output voltage can be obtained:

[0208]

[0209] in,

[0210]

[0211] The corresponding space voltage vector is defined as:

[0212]

[0213] in,

[0214] Since each bridge arm corresponds to three switch states, 27 groups of switch states can be obtained. Substituting them into the defined space voltage vector formula can obtain 27 voltage vectors;

[0215] 4.2: Determination of permanent magnet synchronous generator speed prediction model:

[0216] The dynamic mathematical equation of the permanent magnet synchronous generator is:

[0217]

[0218] in P n is the number of pole pairs of the surface mounted permanent magnet synchronous generator.

[0219] The Laplace transform of the above dynamic mathematical equations gives a first-order model:

[0220]

[0221] in, yes The Laplace transform of is the q-axis current of the surface-mounted permanent magnet synchronous generator at time KT, T m (s) is Tm The Laplace transform of (k), T m (k) is the torque of the surface-mounted permanent magnet synchronous generator at time KT, and s is a complex variable in the frequency domain.

[0222] right Perform Euler discretization to obtain a discrete-time model:

[0223] ω r (k) = ai q (k-1)-bω r (k-1)

[0224] ω r (k+1)=ai q (k)-bω r (k)

[0225] where ω r (k) is the actual speed of the motor at time kT, b=-e -TB / J , the speed prediction model is obtained by rearranging the above formula as follows:

[0226]

[0227] Among them, ω m (k+1) is the predicted output at time (k+1)T, is the given control increment of the q-axis current of the schedule permanent magnet synchronous generator at time kT.

[0228] When the prediction step is extended to N steps, we can get:

[0229]

[0230] Among them, A=-b, B=a; Δi q (k) is the q-axis current control increment of the schedule-mounted permanent magnet synchronous generator at time kT, A i Represents the matrix obtained by raising the matrix A to the i-th power.

[0231] can be simplified to:

[0232] W ωm (k) = W ωr (k)+GΔI q (k)

[0233] in

[0234] Step 5: Establish a discrete permanent magnet synchronous generator current prediction model. The process is as follows:

[0235] According to the coordinate transformation theory, the voltage equation of the surface-mount PMSG mathematical model in the dq coordinate system can be obtained as follows:

[0236]

[0237] Where: u d ,u q The d and q axis voltage components of the stator voltage corresponding to 27 switching states; i d ,i q represents the d and q axis components of the stator current in the dq coordinate system; L s is the stator inductance of the surface-mounted permanent magnet synchronous generator in the dq coordinate system; R s represents the stator resistance; ψ f represents the permanent magnet flux; ω re Represents electrical angular velocity.

[0238] The current state equation is discretized using the forward Euler formula to obtain the discrete permanent magnet synchronous generator current prediction model as follows:

[0239]

[0240] Among them, i d (k),i q (k) represents the stator current component in the dq coordinate system at the current moment; i d (k+1), i q (k+1) is the d,q axis component of the stator current at the next moment; u d (k),u q (k) is the d,q axis component of the optimal voltage vector acting on the system at the last moment in the dq coordinate system; T s is the sampling period;

[0241] Since there is a delay in system control, it is necessary to compensate for the control delay to improve the accuracy of the calculated motor state prediction value to ensure the optimal voltage vector:

[0242]

[0243] Among them, i d (k+1), i q (k+1) is the d,q axis component of the stator current in the dq coordinate system after delay compensation; i d (k+2), i q (k+2) are the predicted values ​​of the d and q axis stator currents at time k+1 respectively; ω re (k+1) is approximately equal to ω re (k).

[0244] Taking into account the uncertainty of system parameters and the influence of external disturbances, the above current equation discretized by Euler is transformed into a vector form:

[0245] i dq (k+1)=Ai dq (k)+BPu(k)+C+B d d(k)

[0246] Among them, among them,

[0247] i dq =[i d i q ] T ,u(k)=[S a S b S c ],

[0248] d(k)=[d d (k) d q (k)] T ,

[0249] Where S a ,S b ,S c Respectively represent the states of the upper and lower switches of the three bridge arms of the inverter, u dc Represents the voltage across the DC bus capacitor.

[0250] When the prediction step is extended to N steps, it is assumed that the rotor speed of the motor remains unchanged within the prediction range, and it is defined as:

[0251]

[0252]

[0253] Where I(k) represents the matrix of the predicted dq-axis current values ​​for the next N time steps at time KT, I * (k) represents the dq-axis reference current matrix for the next N time steps at time KT, A matrix representing the lumped perturbation estimate at time KT for N time steps into the future.

[0254] Thus, we can get:

[0255]

[0256] in

[0257] Step 6, construct the cost function;

[0258] A cost function is constructed to better represent the tracking effect of the predicted value on the expected value. The cost function J1 is as follows:

[0259]

[0260] in, Indicates the reference value of the stator current d, q axis components; i dq (l+1) are the predicted values ​​of d-axis and q-axis stator currents at time l+1, u is the predicted switching frequency, and λ represents the weight coefficient.

[0261] In order to comply with the calculation rules of the sphere decoding algorithm, it is transformed into the following J1 form:

[0262]

[0263] in,

[0264] Bundle Substituting into the cost function J1, we can get:

[0265]

[0266] By further sorting J1 using the matching method, we can get:

[0267] J1=U(k) T HU(k)+2Θ(k) T U(k)+θ(k)=(U(k)+H -1 Θ(k) T H(U(k)+H -1 Θ(k))+const(k) where, const(k) is a constant that is independent of U(k), so the constant term is independent of the switch switching sequence. At this time, the optimization problem can be transformed into:

[0268]

[0269] Since the matrix H is a positive definite symmetric matrix, the equation can be solved to obtain U unc (k)=-H -1 Θ(k), we get the new cost function J'1:

[0270] J'1=(U(k)-U unc (k)) T H(U(k)-U unc (k))

[0271] Among them, U unc(k) is the optimal solution without constraints.

[0272] With the matrix H being a positive definite symmetric matrix, it can be concluded that its inverse matrix is ​​also a positive definite matrix; therefore, there exists a unique reversible lower triangular matrix V that satisfies:

[0273] H -1 =VV T

[0274] make Then the cost function J'1 is transformed into the following formula:

[0275]

[0276] Among them, V -1 is a lower triangular matrix.

[0277] A cost function is constructed to better represent the tracking effect of the predicted value on the expected value. The cost function J2 is as follows:

[0278]

[0279] It is necessary to predict the error in the future optimization time domain and compensate for it.

[0280] e(k)=ω r (k)-ω m (k) and e(k)=e(k+1)=…=e(k+n)

[0281] In order to smoothly transition to the set point, a first-order exponential change is adopted, and the reference trajectory equation is as follows:

[0282] y p (k+j|k)=ω ref (k+j)-α j [ω ref (k)-ω r (k)],j=1…N

[0283] Among them, ω ref (k) is the setting value, α is the weight coefficient, y p (k+1) is the expected output determined by the reference trajectory.

[0284] After prediction correction

[0285] ω' r (k+j|k)=y p (k+j|k)+e(k), j=1…N

[0286] Among them, ω' r (k+j|k) is the corrected reference speed, and e(k) represents the prediction error.

[0287] The above cost function is simplified to:

[0288]

[0289] in,

[0290] make Solving the optimal control rate of the system can be obtained:

[0291] ΔI q (k) = (R + GQG) -1 GQ(W' r -W ωm )

[0292] Among them, Q and R are weight coefficients.

[0293] Step 7, selecting the optimal voltage vector;

[0294] First, the output voltage vector of the three-level inverter is determined by the switching state of the three bridge arms; then, the prediction value at the current moment can be obtained under the action of the prediction model; finally, the optimal voltage vector u is selected according to the designed cost function. opt_1 :

[0295] u opt_1 =argmin J'1

[0296] First, the reference speed at the current moment is collected; then, under the action of the prediction model, the speed prediction value at the current moment can be obtained; finally, the optimal voltage vector u is selected according to the designed cost function. opt_2 :

[0297] u opt_2 =argmin J2

[0298] Step 8, the dq axis differential observer process is as follows:

[0299] Taking into account the uncertainty of system parameters and the influence of external disturbances, the permanent magnet synchronous generator equation in step 5 can be sorted out as follows:

[0300] Let x1(k)=i d (k), then

[0301]

[0302] Among them, d d (k) represents the lumped disturbance, f1(k), f2(k), and f3(k) represent d d (k) first derivative, second derivative, third derivative.

[0303] make Design the corresponding observer, then

[0304]

[0305] Select appropriate observer gains β1, β2, β3, and β4 to achieve asymptotic convergence of the designed dq-axis differential observer. At this time, it is necessary to extract the disturbance information involved in each step of the multi-step FCS-MPC from the observer's estimated value. The specific extraction formula is as follows:

[0306]

[0307] Similarly, a similar observer is designed for the q-axis to obtain the q-axis disturbance information involved in each step of the multi-step FCS-MPC. Based on the disturbance information extracted from the dq-axis difference observer, the required spherical decoding prediction model can be easily obtained.

[0308] Step 9: Introduce state variable d ωl , determine the new state space model;

[0309] In the actual speed regulation process of the motor, the steady-state performance of the speed outer loop will be affected by external disturbances and uncertainty of model parameters. The mechanical motion equation is:

[0310]

[0311] Where, J represents the rotor moment of inertia; T m represents the driving torque; B represents the friction coefficient; ω represents the mechanical angular velocity; T e Represents electromagnetic torque and satisfies

[0312] Due to the uncertainty of model parameters and the influence of external disturbances, the variable d is introduced ωl , transforming the above mechanical motion equation into this form:

[0313]

[0314] in, represents the lumped disturbance on the generator side; It represents the control quantity of the machine-side controller after disturbance compensation; ω represents the mechanical angular velocity.

[0315] Let x 11 =ω,x 12 =d ωl is the new state variable, then the new state space model is expressed as:

[0316]

[0317] Among them, x11 =ω represents the mechanical angular velocity, x 12 =d ωl represents the lumped disturbance, x 13 For x 12 The first derivative of , h1 represents d ωl The second-order differential of It represents the control quantity of the machine-side controller after disturbance compensation;

[0318] Step 10, differential GPIO observer design, the process is as follows:

[0319] According to the new state space model in step 9, the differential GPIO observer is designed in the form of:

[0320]

[0321] in, represents the error variable, represents the estimated value of ω, represents the estimated value of the lumped disturbance, is the estimated value of the first-order derivative of the lumped disturbance; It represents the control quantity of the machine-side controller after disturbance compensation;

[0322] β 11 ,β 12 ,β 13 ,β 14 represents the gain of the extended state observer;

[0323] make The forward Euler formula is used to discretize the above formula to obtain:

[0324]

[0325] At this point, we need to extract β from the estimated value of the differential GPIO observer 11 ,β 12 ,β 13 ,β 14 The perturbation information involved in each step of multi-step speed prediction. The specific extraction formula is as follows:

[0326]

[0327] Step 11: Establish a DC link mathematical model

[0328] The current at the DC side capacitor nodes P, O, and N is expressed as:

[0329]

[0330] i c1 =i pm -ipg

[0331] i c1 +i om =i c2 +i og

[0332] i c2 +i nm =i ng

[0333] Where, C1, C2 represent DC filter capacitors; u c1 ,u c2 Represents the voltage on the DC bus capacitor; i c1 ,i c2 Represents the current flowing through the DC filter capacitor; i pm ,i om ,i nm Indicates the current flowing through nodes P, O, and N on the machine side, i pg ,i og ,i ng Indicates the current flowing to the grid-side nodes P, O, and N.

[0334] Step 12: Establish a network-side mathematical model. The process is as follows:

[0335] The grid-side mathematical model in the dq coordinate system is:

[0336]

[0337] Among them, u d ,u q is the component of the output voltage of the three-level inverter in the d,q coordinate system; e d ,e q is the grid-side voltage component in the d,q coordinate system; i d ,i q is the grid-side current component in the d,q coordinate system; L represents the grid-side filter inductance; R represents the output equivalent resistance; ω ge Represents the grid angular velocity.

[0338] Step 13, sampling and coordinate transformation of grid-side current and voltage;

[0339] In order to achieve effective control of the grid-connected inverter on the grid side and simplify the control system design, it is necessary to convert the information collected in real time in the three-phase stationary coordinate system into the dq coordinate system.

[0340] Clark Transformation:

[0341]

[0342] Park Transform:

[0343]

[0344] Among them, θ ge is the spatial angle of the power grid.

[0345] Step 14, establishing a discrete inner loop power prediction model;

[0346] First, the relationship between the switch state of the three-level inverter and the output voltage vector can be obtained from step 4, and then the collected current, voltage information, and inverter parameter information can be substituted into the prediction model for prediction.

[0347] The grid-side inverter of the wind power generation system is often controlled by grid voltage oriented control. Therefore, the grid-side inverter current equation based on grid voltage vector orientation can be expressed as:

[0348]

[0349] Among them, e d The d-axis voltage representing the grid voltage; ω g Indicates the grid angular frequency.

[0350] According to the grid voltage oriented control and instantaneous power principle, the active power P and reactive power Q of the grid-side inverter can be expressed as:

[0351]

[0352] Among them, e d ,e q ,i d ,i q are the components of the grid voltage and current on the d and q axes respectively.

[0353] Since the grid-side inner loop adopts multi-step power model predictive control, the above power calculation formula needs to be discretized.

[0354] In kT s At any time you can get:

[0355]

[0356] Among them, e d (k),e q (k),i d (k),i q (k) are the d-axis and q-axis components of the grid voltage and current at the current moment, respectively;

[0357] P(k), Q(k) represent the active power and reactive power at the current moment.

[0358] At (k+1)T s At any time you can get:

[0359]

[0360] Among them, e d (k+1),e q (k+1),i d (k+1),i q (k+1) are the d-axis and q-axis components of the predicted values ​​of the grid voltage and current at the next moment respectively; P(k+1) and Q(k+1) represent the predicted values ​​of active power and reactive power at moment k+1 respectively.

[0361] When the sampling time T s Less than When e d (k+1)=e d (k), e q (k+1)=e q (k), then

[0362]

[0363] Therefore

[0364]

[0365] According to the forward Euler formula, the above current equation is discretized to obtain:

[0366]

[0367] Among them, u d (k),u q (k) represents the d and q axis components of the inverter output voltage in the dq coordinate system corresponding to 27 switching states; e d (k) represents the d-axis voltage component of the grid voltage at the current moment;

[0368] The above formula can be rearranged to obtain the multi-step power prediction model:

[0369]

[0370] Among them, u d (k),u q (k) represents the d and q axis components of the inverter output voltage in the dq coordinate system corresponding to the eight switching states.

[0371] Step 15, constructing cost function J3;

[0372] The control goal of the grid side is accurate power tracking and DC side voltage balance, so the cost function J3 can be designed as:

[0373]

[0374] Among them, P * , Q * It represents the reference value of active power and reactive power on the grid side; P(k+1) and Q(k+1) are the predicted values ​​of active power and reactive power at time k+1 respectively.

[0375] Step 16, selecting the optimal voltage vector;

[0376] The predicted value output under each voltage vector is calculated by enumeration, and it is brought into the constructed cost function J3, from which the voltage vector u that minimizes the cost function value is selected. opt_3 ;

[0377]

[0378] Among them, P * , Q * Represents the reference value of active power and reactive power on the grid side; P(k+1) and Q(k+1) are the predicted values ​​of active power and reactive power at time k+1 respectively:

[0379] Find the optimal switching sequence u opt_3 , find the optimal switching sequence within a limited control set.

[0380] u opt_3 =argmin J3

[0381] Step 17, introduce state variable d ul , determine the new state space model;

[0382] Without considering the converter loss, the output power P of the machine-side rectifier is m It can be expressed as:

[0383] P m =u dc i m

[0384] Among them, u dc Represents the DC bus voltage, which can be expressed as u dc =u c1 +u c2 ;i m Indicates the current output from the machine-side converter to the DC bus.

[0385] The current flowing through the DC capacitor is:

[0386]

[0387] Where C represents the DC side capacitance, which can be expressed as C = C1 = C2; i g Represents the current input to the grid-side inverter. The active power P input from the DC side to the grid-side inverter is:

[0388] P=u dc i g

[0389] From the above equation we can get:

[0390]

[0391] Equivalent to

[0392]

[0393] in, represents the total disturbance of the grid-side voltage loop; b u0 It's about b u The estimated value of P * Represents the active power reference value. Let z1=u dc ; z2=d ul , then the new state space model is:

[0394]

[0395] Among them, h2 represents d ul The differential of u0 It's about b u The estimated value of P * Indicates the active power reference value.

[0396] Step 18, designing an extended state observer;

[0397] The extended state observer is designed according to the new state space model in step 17. The grid-side extended state observer is the same as the generator-side extended state observer with an added integral link. The specific form is:

[0398] The improved extended state observer form is as follows:

[0399]

[0400] in, Indicates u dc An estimated value of represents the lumped disturbance d ul The estimated value of ; l1, l2, l3 represent the gains of the improved extended state observer.

[0401] when That is, when the lumped disturbance of the outer loop on the machine side is a constant, and the coefficient matrix of the error state space model is a Hurwitz matrix, the estimated error converges asymptotically to 0, that is, the estimated value can track the actual state asymptotically without error.

[0402] when That is, when the lumped disturbance of the outer loop on the grid side satisfies the time-varying disturbance of the type a1+a2t, the designed extended state observer can achieve error-free asymptotic convergence.

[0403] Step 19, grid-side outer loop control law design, the process is as follows:

[0404] By selecting appropriate observer gains l1, l2, l3, the estimated value of the DC bus voltage can be obtained by the extended state observer designed in step 18. and an estimate of the outer loop lumped disturbance The estimated value obtained by the extended state state observer can be used for controller design, and the specific form is as follows:

[0405]

[0406] in, Indicates the reference value of the voltage outer loop; u u0 Indicates the grid-side controller output; k up represents the controller gain.

[0407] Finally, the algorithm is implemented in Matlab-simulink software, and the simulation results are Figure 5-Figure 11 .

[0408] When the wind speed is changed under the premise of ensuring maximum power tracking, the reference speed of the permanent magnet synchronous generator will change accordingly. Figure 5 , Figure 6 As shown in Figure 2, when the wind speed increases at 0.6 seconds, the speed and electromagnetic torque of the permanent magnet synchronous generator quickly reach the reference value; Figure 7 As shown in , when the wind speed increases at 0.6 seconds, the q-axis current of the permanent magnet synchronous generator can quickly track the reference value to reach a new steady state; Figure 8 As shown in Figure 1, when the wind speed increases at 0.6 seconds, the DC bus voltage is disturbed, but it returns to a steady state in a very short time; Fig. 9 As shown in , when the wind speed increases at 0.6 seconds, the active power on the grid side can well track the active power reference value obtained by the outer loop; Fig.10 As shown in Figure 2, when the grid voltage changes in 0.8 seconds, the DC bus voltage can quickly return to a steady state; Fig.11 As shown in Figure 1, when the disturbance is applied at 0.9 seconds, the speed can quickly return to the reference value.

[0409] The simulation results show that external interference can be well suppressed by the speed outer loop and current inner loop designed above, and the speed can quickly reach the reference value. The voltage loop controller designed above can restore the wind speed and grid voltage to the set value under the influence of disturbance. The dynamic performance and anti-interference performance of the system have been improved to a certain extent.

Claims

1. A multi-step predictive control method for a three-level power generation system based on rolling horizon estimation, characterized in that: The method comprises the following steps: Step 1: Determine the given speed value ω of the machine side speed loop by capturing the wind speed ref ; Step 2: Establish the machine side mathematical model; Step 3: Sample the generator side current and speed and perform Clark-Park transformation. The process is as follows: First, the three-phase current i of the machine-side converter is collected in real time a ,i b ,i c , and then the collected i a ,i b ,i c Perform Clark transformation to obtain the component value i in the α-β coordinate system α ,i β , and perform Park transformation in the α-β coordinate system, the component value i in the α-β coordinate system can be α ,i β Convert to component value i in dq coordinate system d ,i q ; Step 4: Establish a discrete permanent magnet synchronous generator speed prediction model. The process is as follows: 4.1: Determination of output voltage vector of three-level inverter: Assume the three-phase sinusoidal voltage expression is: The inverter output voltage is defined as: but Also U aN +U bN +U cN = 0, then Among them, S a , S b , S c Corresponding to the three switching states of the three bridge arms of the inverter, U aN , U bN , U cN They correspond to the voltage between the inverter output phase A, phase B, phase C and the neutral point N, U oN Represents the inverter output voltage, U dc Represents the voltage across the DC bus capacitor; The relationship between the switching states of the three bridge arms of the three-level inverter and the inverter output voltage is obtained: in, The corresponding space voltage vector is defined as: in, Since each bridge arm corresponds to three switch states, 27 groups of switch states are obtained, which are substituted into the defined space voltage vector formula to obtain 27 voltage vectors; 4.2: Determination of permanent magnet synchronous generator speed prediction model: The dynamic mathematical equation of the permanent magnet synchronous generator is: in P n is the number of pole pairs of the surface-mounted permanent magnet synchronous generator, ψ f is the flux linkage of the surface mounted permanent magnet synchronous generator, T m is the torque of the surface-mounted permanent magnet synchronous generator, J is the moment of inertia, and B is the motor friction coefficient; Laplace transform of the dynamic mathematical equations yields a first-order model: in, yes The Laplace transform of is the q-axis current of the surface-mounted permanent magnet synchronous generator at time KT, T m (s) is T m The Laplace transform of (k), T m (k) is the torque of the surface-mounted permanent magnet synchronous generator at time KT, s is a complex variable in the frequency domain; right Perform Euler discretization to obtain a discrete-time model: ω r (k)=ai q (k-1)-bω r (k-1) ω r (k+1)=ai q (k)-bω r (k) where ω r (k) is the actual speed of the motor at time kT, b=-e -TB / J , the speed prediction model is obtained by rearranging the above formula as follows: Among them, ω m (k+1) is the predicted output at time (k+1)T, is the given control increment of the q-axis current of the schedule-attached permanent magnet synchronous generator at time kT; When the prediction step is extended to N steps, we get: Among them, A=-b, B=a; Δi q (k) is the q-axis current control increment of the schedule-mounted permanent magnet synchronous generator at time kT, A i represents the matrix obtained by raising the matrix A to the power of i; Simplified to: W ωm (k)=W ωr (k)+GΔI q (k) in Step 5: Establish a discrete permanent magnet synchronous generator current prediction model. The process is as follows: According to the coordinate transformation theory, the mathematical model of the surface-mount PMSG in the dq coordinate system is obtained. Cheng Wei: Where: u d ,u q The d and q axis voltage components of the stator voltage corresponding to 27 switching states; i d ,i q represents the d and q axis components of the stator current in the dq coordinate system; L s is the stator inductance of the surface-mounted permanent magnet synchronous generator in the dq coordinate system; R s represents the stator resistance; ψ f represents the permanent magnet flux; ω re represents the electrical angular velocity; The current state equation is discretized using the forward Euler formula to obtain the discrete permanent magnet synchronous generator current prediction model as follows: Among them, i d (k),i q (k) represents the stator current component in the dq coordinate system at the current moment; i d (k+1), i q (k+1) is the d,q axis component of the stator current at the next moment; u d (k),u q (k) is the d, q axis component of the optimal voltage vector acting on the system at the last moment in the dq coordinate system; L s is the stator inductance of the surface-mounted permanent magnet synchronous generator in the dq coordinate system; T s is the sampling period; Since there is a delay in system control, it is necessary to compensate for the control delay to improve the accuracy of the calculated motor state prediction value to ensure the optimal voltage vector: Among them, i d (k+1), i q (k+1) is the d,q axis component of the stator current in the dq coordinate system after delay compensation; i d (k+2), i q (k+2) are the predicted values ​​of the d and q axis stator currents at time k+1 respectively; L s is the stator inductance of the surface-mounted permanent magnet synchronous generator in the dq coordinate system; R s represents the stator resistance; ψ f represents the permanent magnet flux; ω re (k+1) is approximately equal to ω re (k); T s is the sampling period; Taking into account the uncertainty of system parameters and the influence of external disturbances, the above current equation discretized by Euler is transformed into a vector form: i dq (k+1)=Ai dq (k)+BPu(k)+C+B d d(k) Among them, among them, Where S a ,S b ,S c Respectively represent the states of the upper and lower switches of the three bridge arms of the inverter, u dc Indicates the voltage across the DC bus capacitor; When the prediction step is extended to N steps, it is assumed that the rotor speed of the motor remains unchanged within the prediction range, and it is defined as: Where I(k) represents the matrix of the predicted dq-axis current values ​​for the next N time steps at time KT, I * (k) represents the dq-axis reference current matrix for the next N time steps at time KT, The matrix representing the estimated value of the aggregate disturbance at time KT for the next N time steps; Thus, we get: in Step 6, construct the cost function; A cost function is constructed to express the tracking effect of the predicted value on the expected value. The cost function J1 is as follows: in, Indicates the reference value of the stator current d, q axis components; i dq (l+1) are the predicted values ​​of the d and q axis stator currents at time l+1, u is the predicted switching frequency, and λ represents the weight coefficient; In order to comply with the calculation rules of the sphere decoding algorithm, it is transformed into the following J1 form: in, Bundle Substituting into the cost function J1, we get: Using the matching method to further sort out J1, we get: J1=U(k) T HU(k)+2Θ(k) T U(k)+θ(k)=(U(k)+H -1 Θ(k) T H(U(k)+H -1 Θ(k))+const(k) where, const(k) is a constant that is independent of U(k), so the constant term is independent of the switch switching sequence. At this time, the optimization problem is transformed into: Since the matrix H is a positive definite symmetric matrix, the equation can be solved to obtain U unc (k)=-H -1 Θ(k), we get the new cost function J′1: J'1=(U(k)-U unc (k)) T H(U(k)-U unc (k)) Among them, U unc (k) is the optimal solution without constraints; Since the matrix H is a positive definite symmetric matrix, its inverse matrix is ​​also a positive definite matrix; therefore, there exists a unique reversible lower triangular matrix V that satisfies: H -1 =VV T make Then the cost function J′1 is transformed into the following formula: Among them, V -1 is a lower triangular matrix; A cost function is constructed to express the tracking effect of the predicted value on the expected value. The cost function J2 is as follows: Predicting the error in the future optimization time domain and compensating for it can yield: e(k)=ωr(k)-ωm(k) and e(k)=e(k+1)=…=e(k+n) In order to smoothly transition to the set point, a first-order exponential change is adopted, and the reference trajectory equation is as follows: y p (k+j|k)=ω ref (k+j)-a j [oh ref (k)-ω r (k)],j=1…N Among them, ω ref (k) is the set value, α is the weight coefficient, and yp(k+1) is the expected output determined by the reference trajectory; After prediction correction: ω′ r (k+j|k)=y p (k+j|k)+e(k)=j=1…N Among them, ω′ r (k+j|k) is the corrected reference speed, and e(k) represents the prediction error; The above cost function is simplified to: in, make Solving the optimal control rate of the system can be obtained: ΔI q (k)=(R+GQG) -1 GQ(W r '-IN ωm ) Among them, Q, R are weight coefficients; Step 7, selecting the optimal voltage vector; First, the output voltage vector of the three-level inverter is determined by the switching states of the three bridge arms; then, the prediction value at the current moment is obtained under the action of the prediction model; finally, the optimal voltage vector u is selected according to the designed cost function. opt_1 : you opt_1 =angry my J'1 First, the reference speed at the current moment is collected; then, the speed prediction value at the current moment is obtained under the action of the prediction model; finally, the optimal voltage vector u is selected according to the designed cost function. opt_2 : u opt_2 =argmin J2 Step 8, the dq axis differential observer process is as follows: Taking into account the uncertainty of system parameters and the influence of external disturbances, the permanent magnet synchronous generator equation in step 5 is organized as follows: Let x1(k)=i d (k), then Among them, d d (k) represents the lumped disturbance, f1(k), f2(k), and f3(k) represent d d (k) first, second and third derivatives; make Design the corresponding observer, then The observer gains β1, β2, β3, and β4 are selected to make the designed dq-axis differential observer achieve asymptotic convergence; at this time, the disturbance information involved in each step of the multi-step FCS-MPC needs to be extracted from the estimated value of the observer; the specific extraction formula is as follows: Similarly, an observer is designed for the q-axis to obtain the q-axis disturbance information involved in each step of the multi-step FCS-MPC; based on the extracted disturbance information of the dq-axis difference observer, the required spherical decoding prediction model is obtained. Step 9: Introduce state variable d ωl , determine the new state space model; In the actual speed regulation process of the motor, the steady-state performance of the speed outer loop will be affected by external disturbances and uncertainty of model parameters. The mechanical motion equation is: Where, J represents the rotor moment of inertia; T m represents the driving torque; B represents the friction coefficient; ω represents the mechanical angular velocity; T e Represents electromagnetic torque and satisfies Due to the uncertainty of model parameters and the influence of external disturbances, the variable d is introduced ωl , transforming the above mechanical motion equation into this form: in, represents the lumped disturbance on the generator side; represents the control quantity of the machine-side controller after disturbance compensation; ω represents the mechanical angular velocity; Let x 11 =ω,x 12 =d ωl is the new state variable, then the new state space model is expressed as: Among them, x 11 =ω represents the mechanical angular velocity, x 12 =d ωl represents the lumped disturbance, x 13 For x 12 The first derivative of , h1 represents d ωl The second-order differential of It represents the control quantity of the machine-side controller after disturbance compensation; Step 10, differential GPIO observer design, the process is as follows: According to the new state space model in step 9, the differential GPIO observer is designed in the form of: in, represents the error variable, represents the estimated value of ω, represents the estimated value of the lumped disturbance, is the estimated value of the first-order derivative of the lumped disturbance; Represents the control quantity of the machine-side controller after disturbance compensation; β 11 ,β 12 ,β 13 ,β 14 represents the gain of the extended state observer; make The forward Euler formula is used to discretize the above formula to obtain: At this point, we need to extract β from the estimated value of the differential GPIO observer 11 ,β 12 ,β 13 ,β 14 The perturbation information involved in each step of multi-step speed prediction; the specific extraction formula is as follows: Step 11, establishing a DC link mathematical model; Step 12, establishing a network-side mathematical model; Step 13, sampling and coordinate transformation of grid-side current and voltage; Step 14, establishing a discrete inner loop power prediction model; Step 15, constructing a cost function; Step 16, selecting the optimal voltage vector; Step 17, introduce state variable d ul , determine the new state space model; Step 18, designing an extended state observer; Step 19: Design of the grid-side outer loop control law.

2. The multi-step predictive control method for a three-level power generation system based on rolling horizon estimation according to claim 1, characterized in that: In step 2, the specific process is as follows: The mathematical model of the permanent magnet synchronous generator in the dq coordinate system is expressed as: The voltage equation is: Where: u d ,u q represents the d and q axis components of the stator voltage; i d ,i q Indicates the d and q axis components of the stator current; L s is the stator inductance in the dq coordinate system of the surface-mounted permanent magnet synchronous generator, satisfying L s =L d =L q ; R s represents the stator resistance; ω re represents the electrical angular velocity; ψ f represents the permanent magnet flux; The electromagnetic torque equation is: Among them, p n represents the number of pole pairs; T e represents electromagnetic torque; The mechanical motion equation is: Where, ω represents the mechanical angular velocity; J represents the moment of inertia; B represents the friction coefficient; T m Indicates the driving torque.

3. The multi-step predictive control method for a three-level power generation system based on rolling horizon estimation according to claim 1, characterized in that: In step 11, the current at the DC side capacitor nodes P, O, and N is expressed as: i c1 =i pm -i pg i c1 +i om =i c2 +i og i c2 +i nm =i ng Where, C1, C2 represent DC filter capacitors; u c1 ,u c2 Represents the voltage on the DC bus capacitor; i c1 ,i c2 Represents the current flowing through the DC filter capacitor; i pm ,i om ,i nm Indicates the current flowing through nodes P, O, and N on the machine side, i pg ,i og ,i ng Indicates the current flowing to the grid-side nodes P, O, and N.

4. The multi-step predictive control method for a three-level power generation system based on rolling horizon estimation according to claim 1, characterized in that: In step 12, the network-side mathematical model in the dq coordinate system is: Among them, u d ,u q is the component of the output voltage of the three-level inverter in the dq coordinate system; e d ,e q is the grid-side voltage component in the dq coordinate system; i d ,i q is the grid-side current component in the dq coordinate system; L represents the grid-side filter inductance; R represents the output equivalent resistance; ω ge Represents the grid angular velocity.

5. The multi-step predictive control method for a three-level power generation system based on rolling horizon estimation according to claim 1, characterized in that: In step 14, the grid-side inverter adopts a voltage-oriented control method, so the grid-side inverter current equation based on grid voltage vector orientation is expressed as: Among them, u d ,u q is the component of the output voltage of the three-level inverter in the dq coordinate system; e d is the d-axis component of the grid-side voltage; i d ,i q The grid-side current component in the dq coordinate system; L represents the grid-side filter inductance; R represents the output equivalent resistance; ω ge Represents the grid angular velocity.

6. The multi-step predictive control method for a three-level power generation system based on rolling horizon estimation according to claim 1, characterized in that: In step 15, a cost function is constructed, and the cost function J3 is in the following form: Among them, P * ,Q * Indicates the reference value of active power and reactive power; P(k+1), Q(k+1) are (k+1)T s Predicted values ​​of active power and reactive power at each moment.

7. The multi-step predictive control method for a three-level power generation system based on rolling horizon estimation according to claim 6, characterized in that: In step 16, a voltage vector u that minimizes the cost function J2 is selected from the 27 voltage vectors output by the grid-side inverter. opt_3 : u opt_3 =arg min J3。 8. The multi-step predictive control method for a three-level power generation system based on rolling horizon estimation according to claim 1, characterized in that: In step 17, the state variable d is introduced. ul , construct a new state space model; Without considering the converter loss, the output power P of the machine-side rectifier is m It is expressed as: P m =in dc and m Among them, u dc Represents the DC bus voltage, which can be expressed as u dc =u c1 +u c2 ;i m Indicates the current output from the machine-side converter to the DC bus; The current flowing through the DC capacitor is: Wherein, C represents the DC side capacitance, which can be expressed as C=C1=C2; i g Indicates the current input to the grid-side inverter; The active power P input from the DC side to the grid-side inverter is: P=u dc and g From the above equation we get: is equivalent to: in, represents the total disturbance of the grid-side voltage loop; b u0 It's about b u The estimated value of P * Indicates the active power reference value; Let z1=u dc ; z2=d ul , then the new state space model is: Among them, h2 represents d ul The differential of u0 It's about b u The estimated value of P * Indicates the active power reference value.

9. The multi-step predictive control method for a three-level power generation system based on rolling horizon estimation according to claim 1, characterized in that: In step 18, the improved extended state observer of the outer loop is expressed as follows: in, Indicates u dc An estimated value of represents the lumped disturbance d ul The estimated value of; l1, l2, l3 represent the gains of the improved extended state observer; when When the lumped disturbance of the outer loop on the machine side is a constant, and the coefficient matrix of the error state space model is a Hurwitz matrix, the estimated error converges asymptotically to 0, that is, the estimated value can track the actual state asymptotically without error; when That is, when the lumped disturbance of the outer loop on the grid side satisfies the time-varying disturbance of the type a1+a2t, the designed extended state observer can achieve error-free asymptotic convergence.

10. The multi-step predictive control method for a three-level power generation system based on rolling horizon estimation according to claim 1, characterized in that: In step 19, the outer loop control law is designed as follows: Select observer gains l1, l2, l3, and obtain the estimated value of the DC bus voltage using the extended state observer designed in step 18 and an estimate of the outer loop lumped disturbance The estimated value obtained by the extended state observer is used for controller design, and the specific form is as follows: in, Indicates the reference value of the voltage outer loop; u u0 Indicates the grid-side controller output; k up represents the controller gain.

Images