The data in the power system are changing all the time (Li, P. et al., 2022). In the new type of the power system with various new energy sources and emerging loads, the frequency deviation and energy storage charging/discharging state are full of randomness and uncertainty. The manager can foresee a certain range, however, not the exact value, which is similar to the physical state in the quantum field. Inspired by the quantum walks, a quantum model prediction control method has been proposed in this study. The QMPC strategy is mainly composed of three parts: information acquisition, information processing, and instruction output. The overall flowchart of QMPC and the relationship diagram of each part are shown in Figure 1 and Figure 2, respectively.

Quantum walks are the correspondence of classical random walks in the quantum world (Li S. et al., 2023; Pan, L. et al., 2024). Quantum walks can enable the walker to traverse all positions at a faster speed under the influence of quantum coherence and entanglement, thereby accelerating the solution of various problems.

In this study, historical frequency deviation Δ f i values were collected and coded using amplitude. Amplitude code equation is shown in Eq. 1 (Miyamoto and Ueda, 2023): Δ f i ⟩ = p i 0 ⟩ + 1 − p i 1 ⟩ , (1) where Δ f i is the quantum state representation of Δ f i , p i is the amplitude, and 0 and 1 are standard computational basis states.

For the frequency deviation data of the time series, quantum walks in discrete-time one-dimensional position space are chosen for evolution. The system state space of quantum walks is shown in Eq. 2 (Melnikov, A. et al., 2023) H = H W ⊗ H C , (2) where H W is the positional state of the walker, composed of position vectors x ⟩ , x ∈ Z ; x is the integer represented by the current position; H C is a coin state consisting of a linear combination of two basis vectors c ⟩ , c = 0 , 1 ; ⊗ is the direct product operation.

The selection and design of coin operators and walking operators can affect the coherence and entanglement of quantum walks. Therefore, coin operators and walking operators play an important role in quantum walks. Each discrete step of the quantum walk can be divided into two parts: flipping a coin and walking. The coin operator C ∧ is a quantum operator used in quantum computing to model the behavior of a classical coin flip. The commonly used coin operator expression for the Hadamard operation is shown in Eq. 3 (Melnikov, A. et al., 2023) C ∧ = 1 2 1 1 1 − 1 . (3)

The walking operator S ∧ is a description of the translational operations performed by the walker according to the coin state it is in and is used to model the movement of quanta through space. The expression form of the walking operator S ∧ is shown in Eq. 4 (Melnikov, A. et al., 2023) S ∧ = 0 0 ⊗ ∑ x x + 1 ⟨ x + 1 1 ⊗ ∑ x x − 1 x , (4) where ⋅ is the conjugate transpose of ⋅ .

Therefore, one-step evolution of a discrete quantum walk can be described by the unitary operator U = S ∧ C ∧ . At the moment t , the unitary operator U acts on the state ψ ( t − 1 to obtain the evolved state ψ t . As shown in Figure 3, at time t − 1 , the spin-up coin state 1 at grid point x , represented by an upward solid arrow in Figure 3, undergoes a coin flip first under the action of the coin operator C ∧ . Subsequently, combined with the operation of the conditional shift operator S ∧ in Equation 4, the particle at grid point x with coin state 1 moves one step to the right to grid point x + 1 , while the particle with coin state 0 moves one step to the left to grid point x − 1 . At this point, a complete step of quantum walk evolution is achieved.

On repeatedly applying U to the initial state, after walking N steps, the end state of the whole system evolves to ψ N = U N ψ 0 .

Quantum decoding can be accomplished by measuring the amplitude during evolution. Since measurements of quantum states are usually statistical descriptions of classical information, and the expectation values of quantum operators are statistical descriptions of quantum states, the decoding process involves measuring the expectation values of certain quantum operators in order to extract classical information from quantum states. Considering the correlation with frequency deviation, this study chooses to measure the Pauli-Z operator. The estimated value of the frequency deviation at step i in the quantum state Δ f E i can be calculated from the amplitude obtained by decoding (Miyamoto and Ueda, 2023), and the formula is shown in Eq. 5: Δ f E i = 2 ψ i σ Z ψ i − 1 , (5) where ψ i is the quantum state after walking i steps, ψ i is the conjugate transpose of ψ i , and σ Z is the Pauli-Z operator.

In this study, the quantum walks can be divided into three parts: encoding, quantum state walk, and decoding. Amplitude coding is a commonly used coding method in quantum computing. Through amplitude coding, the frequency deviation data can be reflected to the amplitude of the quantum state to establish the mapping relationship between the classical world and the quantum world. Next, the unitary operator U acting on the quantum state is obtained by setting up the walking operator S ∧ and the coin operator C ∧ to specify the walking law. Finally, the expectation value of the Pauli-Z operator is computed during the quantum decoding process to extract the classical information from the quantum state. The pseudocode for quantum walks is as follows:

Model predictive control (MPC) is a closed-loop optimization process based on a predictive model, rolling optimization, and feedback correction. Instead of optimizing the entire future time series at once, the model predictive controller reconsiders the best control strategy for a finite period at each time step. Rolling optimization is one of the characteristics of MPC. Over time, the entire optimization window moves forward, continuously updating and adjusting control inputs to adapt to the current state and objectives of the system, ensuring that the controlled system operates more robustly in dynamic environments. The framework of MPC is shown in Figure 4.

In this study, the system is discretized and the formula of system optimization problem based on model predictive control over a period of time k is shown in Eq. 6 (Hu, J. F. et al., 2021) J x k , U * = ∑ i = 0 P − 1 L x k + i | k , u k + i k + V x k + P k , (6) where J ⋅ is the overall performance metric function and the two variables are the state vector at the current time step x k and the sequence of control inputs U * , where U * = u k k ) , u ( k + 1 k , . . . , u k + P − 1 k ; u k + i | k is the predicted value of the system control input at time step k for the future time step k + i ; L ⋅ is time window performance indicator function; V ⋅ is the terminal cost function; P is the number of steps corresponding to the predicted time domain.

Solve the control sequence U * of length P to obtain the optimal control sequence U opt , take the first element of U opt input to the controlled system, and hold the command until the next time step. Finally, the control quantity is compensated according to the target error caused by disturbances and deviations, the future dynamic output behavior is predicted, and feedback correction is performed to improve the robustness and control accuracy of the system.

In the event of an imbalance between the active power output of generator units and the power demand, the power system experiences power fluctuations, leading to frequency deviations (Zhang X. et al., 2023). The power system model chosen in this study mainly includes conventional units, load, and ESS. The frequency characteristics of the system are shown in Figure 5.

The frequency deviation Δ f and power deviation Δ P can be obtained from Figure 5, the equations are shown in Eqs 7, 8: Δ f = Δ P D + H s , (7) Δ P = Δ P GEN + Δ P BES − Δ P D , (8) where Δ f is frequency deviation; Δ P is power deviation; and Δ P GEN , Δ P BES , and Δ P D are the unit output power, energy storage output power, and load disturbance, respectively. D and H are the load-damping factor and the inertia time constant of the grid, respectively.

In this study, a single-area power grid equipped with a battery energy storage unit is selected as the research object. The primary frequency modulation unit, the secondary frequency modulation unit, and the storage battery model are expressed in the form of the transfer function, and the power system frequency modulation model is obtained as shown in Figure 6.

Battery energy storage (BES) is a common form of the ESS; therefore, BES is chosen for this study. To prevent overcharging and discharging of energy storage batteries, the model takes into account the state of charge of the energy storage battery and inputs it into the quantum model predictive controller through a feedback channel. In practice, the energy storage participates in system frequency regulation with loss and delay in the output power; therefore, the converter and filtering link are replaced by the first-order inertia link equivalent.

In Figure 6, s is the Laplace operator; M is the grid deviation factor; K p and K I are the proportional and integral coefficients of the PI controller, respectively; T G is the governor time constant; T CH , T RH , and F HP are the turbine time constant, re-heater time constant, and re-heater gain, respectively; T BES is the charging and discharging time constant of the BES; C RC is the rated capacity of BES; S O C and S O C 0 are the actual state of charge and the initial state of charge of BES, respectively.

To verify the feasibility and effectiveness of the quantum model prediction strategy proposed in this study, the PID and MPC methods, which are commonly used in power systems, are selected as comparison algorithms. The optimal parameters in the operating range are obtained by empirical methods through continuous debugging. A novel power system that includes a high proportion of energy storage new energy stations is established and simulated on the MATLAB/Simulink platform.

The rated capacity of the conventional unit is set to 500 MW. The ESS consists of 10 battery storage units, and each battery storage unit has a power of 1 MW and a capacity of 0.5 MWh. The base frequency of the system is set to 50 Hz. The initial state of charge (SOC) of the ESS is an important parameter (Zhang et al., 2023c), and the initial SOC is set to 0.6 in this study. The simulation time is set to 60 s. The other main relevant parameters are shown in Table 1.

Step perturbation, which is used to model the effects of sudden load changes on the power system in real situations, is a common type of perturbation encountered in power system simulation. By adding a step disturbance at the 2nd second, the frequency deviation curve, SOC curve, and energy storage output curve of QMPC and comparison methods are shown in Figures 7, 8, 9.

As shown in Figure 7, Figure 8, and Figure 9, the load is suddenly increased at 2 s, resulting in a deficit in the active power of the system. As can be seen from Table 2, the maximum frequency deviation obtained by the PID controller is −0.032 Hz, with a regulation time of 4.34 s. The MPC approach yields a maximum frequency deviation of −0.031 Hz and a regulation time of 3.36 s. The QMPC proposed in this study has the best control effect, with the maximum frequency deviation and the regulation time being −0.022 Hz and 1.80 s, respectively; the maximum frequency deviation and the regulation time are at least 29.1% and 46.4% smaller than the previous two, respectively. Therefore, the QMPC can more effectively restrain the increase in frequency deviation under step disturbance and more quickly stabilize the system in the target state with better dynamic characteristics. This further illustrates the ability of QMPC to better adapt to a sudden load in the power system.

Normally distributed random numbers are generated in MATLAB/Simulink as random perturbations added to the simulation model to simulate load joining with randomness and volatility.

The error integral criterion is commonly used as an evaluation index in the control field. Therefore, in this section, integral absolute error (IAE), integral squared error (ISE), integral time multiple absolute errors (ITAE), and integral time multiple square error (ITSE) are chosen as evaluation indexes. The specific formulas are shown in Eqs 9–12: I A E = ∫ 0 T Δ f t d t , (9) I S E = ∫ 0 T Δ f 2 t d t , (10) I T A E = ∫ 0 T t | Δ f t d t , (11) I T S E = ∫ 0 T t Δ f 2 t d t , (12)

After the random perturbation load is added, the frequency deviation curve, SOC curve, and energy storage output curve of QMPC compared to other comparison methods are shown in Figure 10, Figure 11, and Figure 12, respectively. The error integral indexes, mean absolute error (MAE), and standard deviation (SD) are shown in Table 3.

It can be concluded from Figure 10, Figure 11, Figure 12, and Table 4 that under random perturbation, the maximum frequency deviation obtained by PID, MPC, and QMPC is 0.157 Hz, 0.136 Hz, and 0.069 Hz, respectively, and the average frequency deviation of the absolute values obtained is 0.032 Hz, 0.025 Hz, and 0.016 Hz, respectively. The maximum frequency deviation and average frequency deviation of the absolute values in the proposed method are at least 49.26% and 36.00% smaller than those in the comparison method, respectively. As shown in Table 3, the four integration error indexes of the proposed algorithm have minimum values, with reductions of at least 34.92%, 60.78%, 36.96%, and 63.63%, respectively. The mean absolute error and standard deviation of the proposed method are the minimum values of various algorithms. The data above indicate that the proposed method exhibits stronger overall anti-interference ability and stability. It has excellent control performance and robustness under random disturbance. Therefore, the QMPC is more adaptable to complex operating environments than conventional control methods.