Power management optimization method based on intelligent building green energy system

Abstract

The invention relates to the technical field of intelligent building energy management, in particular to a power management optimization method based on an intelligent building green energy system, which is used for solving the problem of energy supply and demand imbalance caused by high time-varying characteristics of photovoltaic output and load demand in an intelligent building. The method comprises the following steps: firstly, establishing models of an intelligent building energy system, including a photovoltaic system, an energy storage system, a controllable load and a power grid interaction model; performing time-frequency domain analysis on photovoltaic output and load demand by using short-time Fourier transform, and decomposing a time sequence into a trend term, a periodic term and a random fluctuation term; a multi-time scale prediction model is constructed based on the decomposition result, and different prediction methods are adopted for different frequency domain features; and a hierarchical model prediction control framework including a day-ahead planning layer, an hour-level economic dispatching layer and a minute-level real-time control layer is further established, and multi-time-scale coordination control is realized.

Description

Technical Field

[0001] The present invention relates to the technical field of intelligent building energy management, and particularly relates to an optimization method for power management based on an intelligent building green energy system. Background Art

[0002] With the increasingly severe global climate change problem, building energy conservation and emission reduction has become an important way to achieve the "dual carbon" goal. According to statistics, building energy consumption accounts for about 40% of the total social energy consumption, and most of it comes from the consumption of fossil fuels. As an important part of the urban energy system, intelligent buildings are expected to achieve low-carbon and intelligent energy consumption by integrating distributed photovoltaic power generation, energy storage devices, and intelligent energy-consuming equipment. However, in the actual operation process, due to the intermittency and volatility of photovoltaic power generation, as well as the randomness and time-variation of building loads, traditional static energy scheduling methods face many challenges:

[0003] 1. Photovoltaic output volatility: Affected by meteorological conditions, the photovoltaic output power has obvious intra-day fluctuations and seasonal variations. Especially in cloudy weather, the light intensity may change sharply within minutes, resulting in large fluctuations in the photovoltaic output power. According to relevant research data, the photovoltaic output power may fluctuate by up to 80% within 10 minutes in alternating sunny and cloudy weather, seriously affecting the stable operation of the system. In addition, photovoltaic systems in different seasons and different geographical locations exhibit different output characteristics, making it difficult to apply a unified prediction and control method.

[0004] 2. Randomness of load demand: Building loads are affected by various factors such as weather conditions, indoor human activities, and equipment operating states, and have strong randomness and uncertainty. For example, air-conditioning loads are affected by factors such as outdoor temperature, indoor number of people, and building envelope characteristics; lighting loads are affected by sunlight conditions and people's work and rest times; office equipment loads are closely related to working hours and people's behavior habits. The random combination of these factors results in complex and variable time-series characteristics of building loads.

[0005] 3. Lag of traditional scheduling methods: Existing energy scheduling methods are mainly based on hourly predictions and usually use static optimization methods to formulate scheduling plans, making it difficult to cope with minute-level energy supply and demand fluctuations. These methods often rely on historical data for statistical analysis and lack the ability to respond promptly to sudden events (such as weather changes, equipment failures, etc.). In addition, traditional methods usually regard prediction and control as two independent links, and the errors in prediction results cannot be effectively compensated during the control process, resulting in a decline in control performance.

[0006] 4. Multi - time - scale characteristics of the energy system: In the building energy system, there are simultaneously dynamic processes of power electronic devices at the second level, energy supply - demand fluctuations at the minute level, and operation cost optimization problems at the hour level. Traditional methods are difficult to coordinate and handle multi - time - scale problems within a unified framework. This leads to inconsistent, and even conflicting, control objectives at all levels of the system, affecting the overall energy utilization efficiency.

[0007] 5. Challenges in the optimal utilization of energy storage systems: As an important link in balancing photovoltaic power generation and load demand, the charge - discharge strategy of the energy storage system directly affects the economy and reliability of the system. However, characteristics such as limited capacity, imperfect charge - discharge efficiency, and limited cycle life of the energy storage system increase the difficulty of optimal scheduling. How to maximize its economic benefits while ensuring the healthy operation of the energy storage system has become an important challenge for the intelligent building energy system.

[0008] Some solutions have been proposed for the above - mentioned problems in the existing technologies, mainly including:

[0009] Traditional prediction methods: For photovoltaic power generation prediction, currently, physical model methods, statistical methods, and artificial intelligence methods are mainly used. Physical model methods rely on detailed meteorological data and photovoltaic module parameters, with complex calculations and insufficient real - time performance; statistical methods such as time - series analysis (ARIMA, etc.) are difficult to capture non - linear characteristics; while artificial intelligence methods (such as neural networks, support vector machines, etc.) have high accuracy, but have high requirements for training data and insufficient interpretability, making it difficult to perform online parameter adjustment. In terms of load prediction, traditional methods also have deficiencies in accuracy and real - time performance, especially limited prediction ability for short - term fluctuations.

[0010] Traditional scheduling and control methods: Currently, the control strategies of intelligent building energy systems mainly include rule - based control strategies, optimization - based open - loop control strategies, etc. Rule - based control strategies are simple to implement but cannot adapt to complex and changeable operating environments; optimization - based open - loop control strategies can obtain theoretical optimal solutions, but lack an effective feedback adjustment mechanism in the face of real - time disturbances. In recent years, some studies have tried to introduce closed - loop control methods such as model predictive control (MPC), but most of them focus on optimization problems at a single time scale and lack comprehensive consideration of the multi - time - scale characteristics of the system.

[0011] Time - frequency domain analysis methods for energy systems: Traditional time - series analysis methods mainly focus on a single perspective of the time domain or frequency domain and are difficult to comprehensively grasp the time - frequency characteristics of energy data. Although time - frequency joint analysis methods such as wavelet analysis have been applied in energy system analysis, most studies remain in the offline analysis stage and fail to effectively integrate time - frequency analysis results into the prediction and control links, resulting in insufficient information utilization.

[0012] Although model predictive control (MPC) technology has received extensive attention in the field of energy system control in recent years, existing research has mainly focused on the hourly economic dispatch level, and there is less research on the minute-level energy fluctuation response. At the same time, existing methods have deficiencies in dealing with the time-frequency characteristics of photovoltaic power output and load demand, making it difficult to effectively extract and utilize the fluctuation characteristics of multiple time scales, resulting in limited prediction accuracy and control performance. In addition, existing methods usually lack an adaptive parameter adjustment mechanism. Once the system parameters are determined, it is difficult to flexibly adjust them with environmental changes, affecting the long-term stability of the system.

[0013] The relevant research status at home and abroad shows that although certain progress has been made in the prediction and control technology of intelligent building energy systems, there are still obvious deficiencies in coping with highly time-varying photovoltaic power output and load demand. It is mainly reflected in: the lack of a coordinated optimization framework for multiple time scales; the prediction method fails to fully consider the time-frequency characteristics; the control strategy has a lag in response; the system adaptability is insufficient, etc. These problems seriously restrict the in-depth application and efficient utilization of renewable energy in the building field.

[0014] Therefore, there is an urgent need for a real-time optimization method for intelligent building green energy systems that can comprehensively consider the multi-time scale characteristics of energy systems and can quickly respond to minute-level energy supply and demand fluctuations, so as to improve the utilization efficiency of renewable energy, reduce the building operation cost, and enhance the reliability and resilience of system operation. The present invention is an innovative solution proposed in response to the above technical needs. Summary of the Invention

[0015] The present invention aims to solve the technical problems such as response lag, insufficient prediction accuracy, and difficulty in coordinated control of multiple time scales existing in the real-time optimization method of intelligent building green energy systems in the prior art when dealing with the high time-variability of photovoltaic power output and load demand. A power management optimization method based on an intelligent building green energy system is proposed. By combining Fourier time-frequency analysis and a multi-level MPC control strategy, minute-level dynamic optimization scheduling of the building energy system is realized, the local consumption rate of renewable energy is improved, the building energy cost is reduced, carbon emissions are reduced, and the energy resilience and autonomous operation ability of intelligent buildings are enhanced.

[0016] To achieve the above object, the present invention provides a power management optimization method based on an intelligent building green energy system, including the following steps:

[0017] 1. System modeling and state estimation: Establish a model of the intelligent building energy system, including a photovoltaic system, an energy storage system, a controllable load, and a grid interaction model, and update the system state in real time through a state estimation algorithm.

[0018] 1.1 Construction of the photovoltaic system model, considering factors such as solar radiation intensity, ambient temperature, and component characteristics, and establish a photovoltaic output power calculation model;

[0019] 1.2 Energy storage system model construction: Considering factors such as charge-discharge efficiency, capacity decay, and state-of-charge constraints, a dynamic model of the energy storage system is established.

[0020] 1.3 Controllable load model construction, including thermoelectric load models (such as air conditioners, water heaters, etc.) and adjustable power load models (such as electric vehicle charging, some lighting, etc.).

[0021] 1.4 Uncontrollable load model construction: Based on historical data and typical daily characteristics, a statistical model of uncontrollable loads is established.

[0022] 1.5 Grid interaction model construction: Considering factors such as electricity price mechanism and grid connection constraints, an economic model for interaction with the grid is established.

[0023] 1.6 System state estimation algorithm design: Using methods such as Kalman filtering, the system state is updated based on real-time measurement data.

[0024] 2. Fourier time-frequency feature extraction: Using the short-time Fourier transform (STFT) to perform time-frequency domain analysis on photovoltaic output and load demand, extracting frequency domain features on multiple time scales, and decomposing the time series into a trend term, a periodic term, and a random fluctuation term.

[0025] 2.1 Time series data preprocessing, including data normalization, outlier detection and processing, missing value imputation, etc.

[0026] 2.2 Short-time Fourier transform parameter optimization: For the different characteristics of photovoltaic output and load demand, the optimal window function type and window length are determined respectively.

[0027] 2.3 Wavelet packet transform implementation: Selecting a suitable wavelet basis function and determining the decomposition level to further decompose the STFT result.

[0028] 2.4 Multi-time scale frequency domain decomposition: Decomposing the time series into a low-frequency trend term, a medium-frequency periodic term, and a high-frequency fluctuation term.

[0029] 2.5 Frequency domain feature extraction: Calculating statistical features such as energy distribution, phase characteristics, and spectral entropy of each frequency band as the input of the prediction model.

[0030] 3. Multi-time scale prediction model: Based on Fourier time-frequency features, a multi-time scale prediction model for photovoltaic output and load demand is constructed, and different prediction methods are used for different frequency domain features.

[0031] 3.1 Low-frequency trend term prediction model design: Using the support vector regression method, with historical data, weather forecasts, and frequency domain features as inputs, to predict future trend changes.

[0032] 3.2 Design of the intermediate-frequency periodic term prediction model, using the seasonal autoregressive integrated moving average model to capture the periodic change characteristics of the data;

[0033] 3.3 Design of the high-frequency fluctuation term prediction model, using the conditional generative adversarial network to generate multi-scenario prediction results to represent the uncertainty of random fluctuations;

[0034] 3.4 Generation of multi-scenario prediction results, based on the multi-scenario prediction of high-frequency fluctuation terms, combined with the deterministic predictions of low-frequency and intermediate-frequency terms, to generate multi-scenario prediction results of the complete power curve;

[0035] 3.5 Design of the online update mechanism of the prediction model, based on real-time measurement data and prediction errors, dynamically adjust the model parameters to improve the long-term prediction accuracy.

[0036] 4. Hierarchical MPC control framework: Establish a hierarchical MPC control framework including the day-ahead planning layer, the hourly economic dispatch layer, and the minute-level real-time control layer to achieve coordinated control of multiple time scales.

[0037] 4.1 Design of the day-ahead planning layer, based on the day-ahead prediction information, optimize the energy storage charge and discharge plan and the controllable load dispatch plan for the next 24 hours to minimize the daily operating cost;

[0038] 4.2 Design of the hourly economic dispatch layer, based on the rolling updated prediction information, optimize the energy storage operation and controllable load dispatch for the next 1 - 4 hours to minimize the deviation cost;

[0039] 4.3 Design of the minute-level real-time control layer, based on short-term prediction and real-time measurement, optimize the energy storage power and controllable load power for the next 5 - 15 minutes to achieve power balance;

[0040] 4.4 Design of the hierarchical coordination mechanism, through means such as target constraint transfer, prediction information sharing, and result feedback, to ensure the consistency of control decisions at each level;

[0041] 4.5 Design of the multi-objective optimization architecture, comprehensively considering multi-dimensional objectives such as economy, reliability, and environmental protection, and designing reasonable objective functions and constraints;

[0042] 4.6 Consideration of safety constraints, introduce power ramp constraints, energy storage state of charge constraints, equipment operation constraints, etc. to ensure the safe and stable operation of the system.

[0043] 5. Minute-level real-time optimization control: At the minute-level real-time control layer, based on the rolling horizon MPC method, combined with real-time measurement data and short-term prediction results, optimize the energy storage charge and discharge power, controllable load dispatch, and grid interaction power to achieve real-time balance of energy supply and demand.

[0044] 5.1 Design of rolling optimization strategy to determine key parameters such as control period, prediction horizon length, and control horizon length;

[0045] 5.2 Construction of objective function, comprehensively considering multiple sub-objectives such as minimization of power deviation, regulation of state of charge, and control smoothness;

[0046] 5.3 Setting of constraint conditions, including constraints of energy storage system, controllable load, grid interaction, and system power balance;

[0047] 5.4 Design of fast solution algorithm, aiming at the high efficiency requirement of real-time control, developing a fast solution method suitable for online optimization;

[0048] 5.5 Emergency handling mechanism of control strategy, designing fault detection and safety protection strategies to ensure the reliable operation of the system under extreme conditions.

[0049] 6. Adaptive parameter adjustment: Through online learning algorithm, adjust model parameters and control parameters in real time according to prediction error and control performance to improve the adaptability and robustness of the system.

[0050] 6.1 Adaptive adjustment mechanism of prediction model parameters, dynamically update model parameters based on prediction error to improve prediction accuracy;

[0051] 6.2 Adaptive adjustment mechanism of MPC control parameters, dynamically adjust weight coefficients and prediction horizon length according to control performance indicators;

[0052] 6.3 Design of performance evaluation index, establishing a scientific and reasonable evaluation system to quantify prediction and control performance;

[0053] 6.4 Design of online learning algorithm, selecting an efficient learning algorithm suitable for real-time applications to ensure the stability of parameter adjustment;

[0054] 6.5 Design of cold start strategy, solving the problem of insufficient data at the initial operation of the system and accelerating the convergence of system parameters.

[0055] In terms of multi-time scale frequency domain decomposition of photovoltaic output and load demand, the present invention uses short-time Fourier transform to transform the time series signal s(t) into the time-frequency domain:

[0056]

[0057] where S(τ, ω) represents the signal intensity at time τ and frequency ω, s(t) represents the original time series signal, w(t - τ) represents the window function centered at τ, and e -jωt represents the complex exponential basis function with frequency ω.

[0058] For different types of data, the present invention uses different window functions for time-frequency analysis. For photovoltaic output data, the Hanning window is selected as the window function:

[0059] w Hanning (t) = 0.5·[1 - cos(2πt / T)], 0 ≤ t ≤ T;

[0060] where T represents the window length. For load demand data, the Hamming window is selected as the window function:

[0061] w Hamming (t) = 0.54 - 0.46·cos(2πt / T), 0 ≤ t ≤ T;

[0062] By selecting window functions of different widths, frequency domain features on different time scales can be extracted. The present invention uses wavelet packet transform to further decompose the STFT results to obtain multi-scale frequency domain features:

[0063]

[0064] where WT j，k (t) represents the wavelet transform coefficient at scale j and position k, ψ j，k (t) represents the wavelet basis function, and s(τ) represents the original signal.

[0065] The present invention selects the Daubechies wavelet (db4) as the basis function, and its expression is:

[0066] ψ j，k (t) = 2 j / 2 ·ψ(2 j t - k);

[0067] where ψ(t) represents the mother wavelet function, j represents the scale parameter, and k represents the translation parameter.

[0068] Based on the frequency domain decomposition results, the photovoltaic output and load demand are decomposed into components on three time scales:

[0069] P(t) = P trend (t) + P cycle (t) + P fluctuation (t);

[0070] where P(t) represents the original power time series, P trend (t) represents the low-frequency trend term, representing the hourly change trend, P cycle (t) represents the intermediate-frequency periodic term, representing the change characteristics at the scale of dozens of minutes, and P fluctuation (t) represents the high-frequency fluctuation term, representing the random fluctuations at the minute level and below.

[0071] For different frequency domain features, the present invention adopts different prediction methods: for the low-frequency trend term, the support vector regression method is adopted; for the medium-frequency periodic term, the seasonal time series model is adopted; for the high-frequency fluctuation term, the conditional generative adversarial network is used to generate multi-scenario prediction results.

[0072] The expression of the support vector regression method is:

[0073]

[0074] Among them, P ternd (t + k) represents the predicted value of the trend term at time t + k, x represents the input feature vector, and x i represents the feature vector of the i-th training sample. K(x i , x) represents the kernel function, α i and represent the Lagrange multipliers, b represents the bias term, and n represents the number of training samples.

[0075] The expression of the seasonal autoregressive integrated moving average model is:

[0076] Φ(B s )φ(B)(1 - B s ) D (1 - B) d P cycle (t) = Θ(B s )θ(B)ε(t);

[0077] Among them, Φ(B s ) and φ(B) respectively represent the seasonal and non-seasonal autoregressive polynomials, Θ(B s ) and θ(B) respectively represent the seasonal and non-seasonal moving average polynomials, B represents the lag operator, s represents the seasonal period, D and d respectively represent the seasonal and non-seasonal differencing orders, and ε(t) represents white noise.

[0078] The training objective of the conditional generative adversarial network is:

[0079]

[0080] Among them, x represents the real sample, z represents the random noise, y represents the conditional variable, p data (x) represents the real data distribution, p z (z) represents the noise distribution, G(z|y) represents the sample generated based on the condition y and the noise z, and D(x|y) represents the discrimination result of the discriminator for the sample x under the condition y.

[0081] The output of the multi-time scale prediction model serves as the input to the hierarchical MPC control framework. In the minute-level real-time control layer, the expression of the MPC optimization problem is:

[0082]

[0083] Constraints:

[0084] x(t+i+1|k) = Ax(t+i|k) + Bu(t+i|k) + Ed(t+i|k);

[0085] y(t+i|k) = Cx(t+i|k);

[0086] u min ≤ u(t+i|k) ≤ u max ;

[0087] x min ≤ (t+i|k) ≤ x max ;

[0088] Δu min ≤ u(t+i|k) - u(t+i - 1|k) ≤ Δu max ;

[0089] where x(t+i|k) represents the state prediction at time t+i at time k, u(t+i|k) represents the control input, d(t+i|k) represents the predicted disturbance, y(t+i|k) represents the system output, N represents the prediction horizon length, Q and R represent the state weighting matrix and the control input weighting matrix respectively, A, B, E, and C are system model matrices, u min 、u max 、x min 、x max 、Δu min and Δu max represent the constraints of the control input, state variables, and their rates of change respectively.

[0090] The adaptive adjustment of the prediction model parameters is based on an online learning algorithm, and the update formula is:

[0091]

[0092] where θ (t) represents the model parameters at time t, η represents the learning rate, represents the gradient of the loss function with respect to the parameters.

[0093] One of the adaptive adjustment methods for MPC control parameters is the update of the weight coefficients Q and R:

[0094] Q (t+1) = Q (t) + α·[(xT x) (t) -(x T x) target ;

[0095] R (t+1) = R (t) + β·[(u T u) (t) -(u T u) target ;

[0096] Wherein, Q (t) and R (t) respectively represent the state weighting matrix and the control input weighting matrix at time t, α and β represent adjustment coefficients, (x T x) (t) and (u T u) (t) respectively represent the sum of squares of states and the sum of squares of control inputs at time t, (x T x) target and (u T u) target represent the target values.

[0097] The present invention has the following beneficial effects:

[0098] 1. Improve prediction accuracy: Through Fourier time-frequency feature extraction, the photovoltaic output and load demand are decomposed into frequency-domain features on different time scales, and different prediction methods are adopted accordingly, significantly improving the prediction accuracy.

[0099] 2. Achieve minute-level real-time optimization: Based on the hierarchical MPC control framework, the present invention realizes multi-time-scale coordinated optimization from day-ahead planning to minute-level real-time control, can quickly respond to minute-level energy supply and demand fluctuations, enables the system to better cope with random disturbances such as cloud changes, sunny and cloudy changes, and load mutations, and improves the dynamic performance of the system. Detailed implementation manners

[0100] The present invention will be further described in detail below.

[0101] 1. System architecture:

[0102] The intelligent building green energy system of the present invention includes a photovoltaic power generation system, an energy storage system, a controllable load, an uncontrollable load, and a grid interaction interface. The core control modules of the system include: a data acquisition and preprocessing module, a Fourier time-frequency feature extraction module, a multi-time-scale prediction module, a hierarchical MPC control module, and an adaptive parameter adjustment module.

[0103] 1.1 Hardware architecture:

[0104] The hardware architecture of the system mainly includes the following parts:

[0105] 1. Sensing and measurement devices: including a photovoltaic system power meter, an energy storage system power meter and state of charge monitoring device, a load power meter, indoor and outdoor environmental sensors (temperature, humidity, light, etc.), a grid interaction power meter, etc.

[0106] 2. Communication network: Adopting a multi-level communication architecture, the underlying devices are connected to the field controller through industrial buses such as RS-485 and Modbus; the field controller is connected to the central control system through Ethernet; the central control system is connected to the cloud platform through the Internet to achieve remote monitoring and data analysis.

[0107] 3. Control and execution devices: including energy storage system power converters (PCS), controllable load controllers (such as variable frequency air conditioner controllers, smart sockets, etc.), relays, switches and other execution devices.

[0108] 4. Computing and storage devices: including edge computing servers (for real-time data processing and control computing), cloud servers (for big data analysis and remote monitoring), data storage devices, etc.

[0109] 1.2 Software architecture:

[0110] The software architecture of the system adopts a hierarchical design and mainly includes the following layers:

[0111] 1. Data acquisition layer: Responsible for collecting real-time data from various sensors, measurement devices and controllers, and realizing data standardization, time synchronization and preliminary screening.

[0112] 2. Data preprocessing layer: Performs operations such as cleaning, anomaly detection, missing value processing, and data compression on the collected raw data to improve data quality.

[0113] 3. Feature extraction layer: Realizes the Fourier time-frequency feature extraction function, and converts the original time series into multi-time scale frequency domain features.

[0114] 4. Prediction and analysis layer: Based on the extracted features, runs a multi-time scale prediction model to generate prediction results of photovoltaic output and load demand.

[0115] 5. Optimization control layer: Based on the prediction results and system status, runs a hierarchical MPC control algorithm to generate an optimization control strategy.

[0116] 6. Execution layer: Converts the control strategy into specific control instructions, issues them to each execution device, and monitors the execution status.

[0117] 7. Visualization and interaction layer: Provides a human-computer interaction interface to realize functions such as system monitoring, parameter setting, and report generation.

[0118] 8. Security management layer: Responsible for the data security, device security, and operation security management of the system, including functions such as access control, data encryption, and fault detection.

[0119] Each layer conducts data exchange and function calls through standardized interfaces to ensure the modularity and scalability of the system.

[0120] 2. Intelligent building energy system model:

[0121] 2.1 Photovoltaic system model:

[0122] The output power of the photovoltaic system can be expressed as:

[0123] P PV (t) = η PV ·S·G(t)·[1 - β(T c (t) - T ref )];

[0124] Among them, P PV (t) represents the photovoltaic output power at time t (kW), η PV represents the conversion efficiency of the photovoltaic module, S represents the total area of the photovoltaic module (m 2 ), G(t) represents the solar radiation intensity at time t (kW / m 2 ), β represents the temperature coefficient (usually 0.004 - 0.006 / °C), T c (t) represents the component temperature at time t (°C), T ref represents the reference temperature under standard test conditions (usually 25°C).

[0125] The component temperature Tc ( t) can be calculated by the following formula:

[0126]

[0127] Among them, T a (t) represents the ambient temperature at time t (°C), and NOCT represents the nominal operating cell temperature (usually 42 - 48°C).

[0128] Considering the dynamic characteristics of the photovoltaic system, a first-order lag link is introduced to describe the dynamic process from the change in radiation intensity to the change in output power:

[0129]

[0130] Among them, P PV,real (t) represents the actual output power (kW), P PV (t) represents the theoretical output power (kW) calculated according to the static model, τ PVRepresents the time constant (usually 10 - 30 seconds).

[0131] In addition, considering the efficiency characteristics of the PV inverter, the actual output power needs to be multiplied by the inverter efficiency:

[0132] P PV，out (t) = η inv (P PV，real (t)) · P PV，real (t);

[0133] where, T PV，out (t) represents the inverter output power (kW), η inv (P PV，real (t)) represents the inverter efficiency, which is a function of power and is usually represented by a piecewise linear function.

[0134] 2.2 Energy storage system model:

[0135] The dynamic model of the energy storage system can be expressed as:

[0136]

[0137] where, SOC(t) represents the state of charge of the energy storage at time t, η ch represents the charging efficiency, η dis represents the discharging efficiency, P ch (t) represents the charging power (kW) at time t, P dis (t) represents the discharging power (kW) at time t, Δt represents the time step (h), and E cap represents the rated capacity (kWh) of the energy storage system.

[0138] The operating constraints of the energy storage system include:

[0139] 0 ≤ P ch (t) ≤ P ch，max ;

[0140] 0 ≤ P dis (t) ≤ P dis，max ;

[0141] P ch (t) · P dis (t) = 0;

[0142] SOC min ≤ SOC(t) ≤ SOC max ;

[0143] where, P ch，max represents the maximum charging power (kW), P dis，max represents the maximum discharging power (kW), SOC min [[ID=8,5]]and SOCmax respectively represent the minimum and maximum allowable state of charge.

[0144] Consider the charge and discharge power ramp rate constraints:

[0145] -r down ≤P ch (t) - P ch (t - 1) ≤ r up ;

[0146] -r down ≤P dis (t) - Pdis(t - 1) ≤ r up ;

[0147] where r down and r up respectively represent the maximum downward and upward rates (kW / min).

[0148] Consider the impact on the life of the energy storage system and introduce a capacity degradation model:

[0149]

[0150] where E cap (t) represents the effective capacity (kWh) at time t, α deg represents the capacity degradation coefficient, DOD(t) represents the depth of discharge, N cycle (t) represents the equivalent cycle number, β temp represents the temperature coefficient, T bat (t) represents the battery temperature (°C), T ref，bat represents the reference temperature (°C).

[0151] 2.3 Load model:

[0152] The building load can be divided into two parts: controllable load and uncontrollable load:

[0153] P load (t) = P ctrl (t) + P non-ctrl (t);

[0154] where P load (t) represents the total load power (kW) at time t, P ctrl (t) represents the controllable load power (kW), P non-ctrl (t) represents the uncontrollable load power (kW).

[0155] The controllable load mainly includes thermoelectric devices such as air conditioners and water heaters, and its model can be expressed as:

[0156]

[0157] Among them, T in (t) represents the indoor temperature (°C) at time t, C th represents the heat capacity of the building (kWh / °C), Q env (t) represents the environmental heat exchange power (kW) at time t, and COP represents the energy efficiency coefficient of the equipment.

[0158] The environmental heat exchange power Q env (t) can be expressed as:

[0159] Q env (t) = UA·(T out (t) - T in (t)) + Q int (t) + Q sol (t);

[0160] Among them, UA represents the overall heat transfer coefficient of the building (kW / °C), T out (t) represents the outdoor temperature (°C), Q int (t) represents the internal heat gain (kW), Q sol (t) represents the solar radiation heat gain (kW).

[0161] The operation constraints of the controllable load include:

[0162] 0 ≤ P ctrl (t) ≤ P ctrl，max ;

[0163] T min ≤ T in (t) ≤ T max ;

[0164] Among them, P ctrl，max represents the maximum controllable load power (kW), T min and T max respectively represent the lower and upper limits of the indoor temperature (°C).

[0165] Considering the thermal comfort index, the predicted mean vote (PMV) is introduced:

[0166] PMV(t) = f(T in (t), RH(t), v air (t), MET, CLO);

[0167] Among them, RH(t) represents the indoor relative humidity, v air (t) represents the air velocity (m / s), MET represents the metabolic rate, and CLO represents the clothing thermal resistance.

[0168] The thermal comfort constraint is:

[0169] -0.5 ≤ PMV(t) ≤ 0.5;

[0170] For non-controllable loads, a statistical model is used for description:

[0171] P non-ctrl (t) = P base (d, h) + P weather (t) + P random (t);

[0172] Among them, P base (d, h) represents the base load based on date d and hour h, and P weather (t) represents the load adjustment amount related to weather, and P random (t) represents the random fluctuation term.

[0173] 2.4 Grid interaction model:

[0174] Assume that the power balance relationship between the building and the grid is:

[0175] P grid (t) = P load (t) - P PV (t) - P dis (t) + P ch (t);

[0176] Among them, P grid (t) represents the electricity purchased from the grid at time t (kW, positive value indicates electricity purchase, negative value indicates electricity sale).

[0177] The constraint conditions for grid interaction are:

[0178] P grid,min ≤ P grid (t) ≤ P grid,max ;

[0179] Among them, P grid,min and P grid，max represent the maximum electricity sale power and the maximum electricity purchase power (kW), respectively.

[0180] Considering the power factor constraint:

[0181] PF min ≤ PF(t) ≤ 1;

[0182] Among them, PF(t) represents the power factor at time t, and PF min represents the minimum power factor requirement (usually 0.9 or 0.95).

[0183] The grid electricity price model includes two parts: the electricity volume price and the demand price:

[0184]

[0185] Among them, C grid (t) represents the grid cost (yuan), C energy (t) represents the electricity price (yuan / kWh), C demand represents the demand price (yuan / kW), T bill represents the billing period.

[0186] The electricity sales revenue model is:

[0187] R sell (t) = C sell (t)·(-P grid (t))·Δt·I(P grid (t) < 0);

[0188] Among them, R sell (t) represents the electricity sales revenue (yuan), C sell (t) represents the electricity sales price (yuan / kWh), I(·) represents the indicator function.

[0189] 3. Fourier time-frequency feature extraction:

[0190] 3.1 Data preprocessing:

[0191] Before performing Fourier time-frequency feature extraction, it is necessary to preprocess the original time series data, which mainly includes the following steps:

[0192] 1. Data cleaning: Detect and process problems such as outliers and missing values to ensure the integrity and validity of the data.

[0193] 2. Data normalization: Unify data with different dimensions to the same scale for subsequent processing. The normalization formula is:

[0194]

[0195] Among them, x norm represents the normalized data, x represents the original data, μ represents the mean, and σ represents the standard deviation.

[0196] 3. Trend removal: For time series with obvious trends, first perform trend removal to better extract periodic fluctuation features. Trend removal can use difference or polynomial fitting methods.

[0197] 4. Data segmentation: Divide the long time series data into several overlapping short time windows for short-time Fourier transform.

[0198] 3.2 Short-time Fourier transform:

[0199] To capture the time-varying characteristics of photovoltaic output and load demand, the present invention uses the short-time Fourier transform for time-frequency analysis. The calculation formula of STFT is:

[0200]

[0201] where S(τ, ω) represents the signal intensity at time τ and frequency ω, s(t) represents the original time series signal, w(t - τ) represents the window function centered at τ, and e -jwt represents the complex exponential basis function with frequency ω.

[0202] In actual calculations, the discrete short-time Fourier transform (DSTFT) is used:

[0203]

[0204] where S(m, k) represents the STFT coefficient at time point m and frequency point k, and N represents the window length.

[0205] In the present invention, appropriate window functions and window lengths are selected respectively according to the different characteristics of photovoltaic output and load demand. For photovoltaic output, the Hanning window is used and a window length of 10 minutes is selected; for load demand, the Hamming window is used and a window length of 15 minutes is selected.

[0206] The expression of the Hanning window is:

[0207] w Hanning (n) = 0.5·[1 - cos(2πn / (N - 1))], 0 ≤ n ≤ N - 1;

[0208] The expression of the Hamming window is:

[0209] w Hamming (n) = 0.54 - 0.46·cos(2πn / (N - 1)), 0 ≤ n ≤ N - 1;

[0210] The amplitude of the time-frequency spectrum can be calculated by the following formula:

[0211]

[0212] The phase spectrum can be calculated by the following formula:

[0213] The phase spectrum can be calculated by the following formula:

[0214]

[0214] 3.3 Wavelet packet transform:

[0215] To further extract the frequency-domain characteristics of multiple time scales, the present invention uses wavelet packet transform to decompose the STFT results:

[0216]

[0217] Among them, WT j，k (t) represents the wavelet transform coefficient at scale j and position k, and ψ j，k (t) represents the wavelet basis function, and s(τ) represents the original signal.

[0218] The discrete form of the wavelet packet transform is as follows:

[0219]

[0220] The present invention selects the Daubechies wavelet (db4) as the basis function and adopts a three-level wavelet packet decomposition. The scale function coefficients of the Daubechies wavelet (db4) are:

[0221]

[0222] The wavelet function coefficients are:

[0223] g0 = h3, g1 = -h2, g2 = h1, g3 = -h0;

[0224] The calculation of the wavelet packet decomposition adopts the filter bank method to perform high-pass and low-pass filtering operations on the signal:

[0225]

[0226] Among them, A j，k represents the low-frequency sub-band coefficient, and D j，k represents the high-frequency sub-band coefficient.

[0227] 3.4 Multi-time scale frequency domain decomposition:

[0228] Based on the results of STFT and wavelet packet transform, the photovoltaic output and load demand are decomposed into components of three time scales:

[0229] P(t) = P trend (t) + P cycle (t) + P fluctuation (t);

[0230] Among them, P(t) represents the original power time series, P trend (t) represents the low-frequency trend term (corresponding to the hourly change trend), P cycle (t) represents the intermediate-frequency periodic term (corresponding to the change characteristics of dozens of minutes), P fluctuation (t) represents the high-frequency fluctuation term (corresponding to the random fluctuations at the minute level and below).

[0231] The specific process of the frequency domain decomposition is as follows:

[0232] (1) Perform STFT on the original time series to obtain the time-frequency spectrum;

[0233] (2) Perform wavelet packet decomposition on the STFT result to obtain the coefficients of different frequency bands;

[0234] (3) According to the frequency band range, divide the wavelet packet coefficients into three parts: low frequency (0 - 0.2 mHz), medium frequency (0.2 - 2 mHz), and high frequency (> 2 mHz);

[0235] (4) Perform inverse transformation on the coefficients of each frequency band to obtain the corresponding time-domain signals P trend (t), P cycle (t), and P fluctuation (t).

[0236] Table 1 lists the characteristic parameters of the three time-scale components:

[0237] Component Frequency range Time scale Main features Energy proportion <![CDATA[P trend (t)]]> 0 - 0.2 mHz Hourly level Daily variation trend, weather influence 65％～75％ <![CDATA[P cycle (t)]]> 0.2 - 2 mHz Ten - minute level Equipment start - stop, periodic behavior 20％～30％ <![CDATA[P fluctuation (t)]]> > 2 mHz Minute level and below Random disturbance, measurement noise 5％～10％ ;

[0238] 3.5 Frequency-domain feature extraction:

[0239] For each frequency-domain component obtained by decomposition, calculate the following characteristic parameters as the input of the prediction model:

[0240] 1. Energy distribution: The energy proportion of each frequency band, and the calculation formula is:

[0241]

[0242] Among them, E i represents the energy proportion of the i-th frequency band, and Ω i represents the frequency range of the i-th frequency band.

[0243] 2. Spectrum entropy: Characterize the uncertainty of the spectrum distribution, and the calculation formula is:

[0244]

[0245] Among them, represents the energy proportion of frequency w.

[0246] 3. Spectrum centroid: Characterize the "center of gravity" position of the spectrum, and the calculation formula is:

[0247]

[0248] 4. Spectrum bandwidth: Characterize the dispersion degree of the spectrum, and the calculation formula is:

[0249]

[0250] 5. Spectrum flatness: Characterize the smoothness of the spectrum, and the calculation formula is:

[0251]

[0252] Among them, N represents the number of frequency points.

[0253] 6. Spectrum attenuation rate: Characterizes the attenuation speed of high-frequency components, and the calculation formula is:

[0254]

[0255] In addition, the cross features in the time-frequency domain are also calculated, such as the energy ratio and phase difference between different frequency bands, to capture the inter-band relationships.

[0256] 4. Multi-time scale prediction model:

[0257] For different frequency domain features, the present invention adopts different prediction methods:

[0258] 4.1 Low-frequency trend term prediction:

[0259] For the low-frequency trend term P trend [[ID=2-eight]](t), a support vector regression method is used for prediction. The expression of support vector regression is:

[0260]

[0261] Among them, P trend (t + k) represents the predicted value of the trend term at time t + k, x represents the input feature vector, and x i represents the feature vector of the i-th training sample, K(x i , x) represents the kernel function, α i and represent the Lagrange multipliers, b represents the bias term, and n represents the number of training samples.

[0262] The present invention selects the radial basis function (RBF) as the kernel function:

[0263] K(x i , x) = exp(-γ||x i [[ID=5-four]] - x|| 2 ;

[0264] Among them, γ is the kernel parameter.

[0265] The optimization problem of support vector regression is expressed as:

[0266]

[0267] Constraint conditions:

[0268] y i - w T φ(xi ) - b ≤ ε + ξ i ;

[0269]

[0270] Among them, w represents the weight vector, φ(x i ) represents the feature mapping, ε represents the width of the insensitive region, ξ i and represent the slack variables, and C represents the penalty parameter.

[0271] The input features include: historical power data, weather forecast information (temperature, radiation intensity, cloud cover, etc.), time features (hour, date type, etc.), and frequency domain features (energy distribution, spectral entropy, etc.). Feature selection adopts the recursive feature elimination (RFE) method to retain the most important feature subset. Hyperparameter optimization adopts the grid search and cross-validation methods.

[0272] For the case with a relatively high degree of non-linearity, an ensemble learning method such as Random Forest or Gradient Boosting Tree can also be adopted to improve the prediction accuracy.

[0273] 4.2 Medium-frequency periodic term prediction:

[0274] For the medium-frequency periodic term P cycle (t), a seasonal autoregressive integrated moving average model is adopted:

[0275] Φ(B s )φ(B)(1 - B<� s ) D (1 - B) d P cycle (t) = Θ(B s )θ(B)ε(t);

[0276] Among them, Φ(B s ) and φ(B) represent the seasonal and non-seasonal autoregressive polynomials respectively, Θ(B s ) and θ(B) represent the seasonal and non-seasonal moving average polynomials respectively, B represents the lag operator, s represents the seasonal period, D and d represent the seasonal and non-seasonal differencing orders respectively, and ε(t) represents white noise.

[0277] Specifically, the expressions of φ(B), θ(B), Φ(B s ) and Θ(B s ) are:

[0278] φ(B) = 1 - φ1B - φ2B 2 - … - φ p B It should be noted that there is a possible error in the original text where "<� s " is used. It might be a misrepresentation. I've translated it as is while pointing it out for your reference.p ;

[0279] θ(B) = 1 + θ1B + θ2B 2 + … + θ q B q ;

[0280] Φ(B s ) = 1 - Φ1B s - Φ2B 2s - … - Φ P B Ps ;

[0281] Θ(B s ) = 1 + Θ1B s + Θ2B 2s + … + Θ Q B Qs ;

[0282] where p, q, P, and Q represent the non-seasonal autoregressive order, non-seasonal moving average order, seasonal autoregressive order, and seasonal moving average order, respectively.

[0283] According to the different periodic characteristics of photovoltaic output and load demand, appropriate SARIMA model parameters are determined respectively. For photovoltaic output, the SARIMA(2, 1, 1)(1, 1, 1)48 model is adopted (assuming the data acquisition interval is 30 minutes, then there are 48 points in a day); for load demand, the SARIMA(3, 1, 2)(1, 1, 1)48 model is adopted.

[0284] The model parameter estimation uses the maximum likelihood estimation (MLE) method, and the optimization objective is:

[0285]

[0286] where L represents the likelihood function.

[0287] The model selection is based on the Akaike information criterion (AIC) and the Bayesian information criterion (BIC):

[0288] AIC = -2ln(L) + 2k;

[0289] BIC = -2ln(L) + k ln(n);

[0290] where k represents the number of model parameters and n represents the number of samples.

[0291] 4.3 High-frequency fluctuation term prediction:

[0292] For the high-frequency fluctuation term P fluctuation(t), a conditional generative adversarial network is used to generate multi-scenario prediction results. The CGAN consists of a generator G and a discriminator D, and its training objective is:

[0293]

[0294] Among them, x represents the real sample, z represents the random noise, y represents the conditional variable, and p data (x) represents the real data distribution, and p z (z) represents the noise distribution. G(z|y) represents the sample generated based on the conditional variable y and the noise z, and D(x|y) represents the discrimination result of the discriminator for the sample x under the condition y.

[0295] The conditional variable y includes historical high-frequency fluctuation data, the current system state, and low-frequency and medium-frequency prediction results. Both the generator and the discriminator adopt the deep convolutional neural network (DCNN) structure.

[0296] The network structure of the generator G is: input layer (noise vector z and conditional variable y) → fully connected layer → reshape layer → transposed convolutional layer 1 (64 5×5 convolutional kernels, stride 2) → batch normalization layer → LeakyReLU activation → transposed convolutional layer 2 (32 5×5 convolutional kernels, stride 2) → batch normalization layer → LeakyReLU activation → transposed convolutional layer 3 (1 5×5 convolutional kernel, stride 1) → tanh activation → output layer.

[0297] The network structure of the discriminator D is: input layer (sample x and conditional variable y) → convolutional layer 1 (32 5×5 convolutional kernels, stride 2) → LeakyReLU activation → convolutional layer 2 (64 5×5 convolutional kernels, stride 2) → batch normalization layer → LeakyReLU activation → convolutional layer 3 (128 5×5 convolutional kernels, stride 2) → batch normalization layer → LeakyReLU activation → fully connected layer → sigmoid activation → output layer.

[0298] The network parameters are optimized using the Adam optimizer, with a learning rate of 0.0002, β1 = 0.5, and β2 = 0.999. The training process includes two stages: pre-training and fine-tuning. Historical data is used for pre-training, and the data in the most recent period is used for fine-tuning.

[0299] After the training is completed, by randomly sampling the noise vector z, multiple possible scenarios of high-frequency fluctuations can be generated, which characterize the uncertainty of the prediction.

[0300] 4.4 Fusion of multi-scenario prediction results:

[0301] The prediction results of the three time scales are fused to obtain the final predicted power:

[0302]

[0303] Among them, represents the predicted value of the power at time t + k, and respectively represent the predicted values of the three time-scale components.

[0304] For high-frequency fluctuation terms, since the CGAN generates multiple scenarios, scenario probability weighting is required:

[0305]

[0306] Among them, S represents the number of scenarios, and π s represents the probability weight of the s-th scenario, represents the predicted value of high-frequency fluctuations under the s-th scenario.

[0307] In addition to the expected value, the probability distribution characteristics of the prediction, such as quantiles, standard deviations, etc., are also calculated to provide a basis for risk assessment and robust control.

[0308] In the post-processing stage of the prediction results, physical constraint conditions are also applied, such as the photovoltaic output power cannot be less than zero and cannot exceed the rated power, etc., to ensure the physical rationality of the prediction results.

[0309] 5. Hierarchical MPC control framework:

[0310] The present invention adopts a hierarchical MPC control framework, which includes three levels: the day-ahead planning level, the hourly economic dispatch level, and the minute-level real-time control level. The control objectives, decision variables, constraint conditions, and time scales of the three levels are shown in Table 2:

[0311] Control level Main objective Decision variable Time scale Rolling period Solution algorithm Day - ahead planning layer Minimize daily operating cost Daily charge - discharge plan of energy storage 24 hours 1 day Mixed - integer linear programming Hourly scheduling layer Minimize deviation cost Hourly adjustment of energy storage, controllable load planning 1 - 4 hours 15 minutes Quadratic programming Minute - level control layer Minimize power deviation Minute - level control of energy storage, real - time load adjustment 10 - 15 minutes 1 minute Model predictive control ;

[0312] 5.1 Day-ahead planning level:

[0313] Based on the day-ahead prediction information, the day-ahead planning level optimizes the energy storage charge and discharge plan for the next 24 hours, and its optimization objective is to minimize the daily operating cost:

[0314]

[0315] The constraint conditions include:

[0316] P grid，buy (t) - P grid，sell (t) = P load，pred (t) - P PV，pred (t) - P dis (t) + P ch (t);

[0317] 0 ≤ P grid，buy (t) ≤ Pgrid，max ;

[0318] 0 ≤ P grid，sell (t) ≤ P rid，max ;

[0319] P grid，buy (t) · P grid，sell (t) = 0;

[0320] 0 ≤ P ch (t) ≤ P ch，max ;

[0321] 0 ≤ P dis (t) ≤ P dis，max ;

[0322] P ch (t) · P dis (t) = 0;

[0323]

[0324] SOC min ≤ SOC(t) ≤ SOC max ;

[0325] SOC(T) = SOC target ;

[0326] where C grid,buy (t) and C grid，sell (t) represent the grid power purchase and power selling prices (yuan / kWh) at time t respectively, P grid,buy (t) and P grid，sell (t) represent the power purchase and power selling powers (kW) at time t respectively, C bat represents the depreciation cost of battery charge and discharge (yuan / kWh), P load，pred (t) and P PV,pred (T) represent the load and PV predicted power (kW) at time t respectively, T represents the planning time domain length (24 hours), SOC targt represents the target state of charge.

[0327] For the consideration of demand-based electricity price, introduce an auxiliary variable P demand representing the maximum power within the billing period:

[0328]

[0329] Modify the objective function to:

[0330]

[0331] Considering the influence of energy storage life, introduce the cycle aging cost:

[0332]

[0333] Among them, α cycle represents the cost per cycle (yuan / cycle), and E cap represents the rated capacity (kWh).

[0334] The updated objective function is as follows:

[0335]

[0336] 5.2-hour economic dispatch layer:

[0337] The hourly economic dispatch layer optimizes the energy storage charge and discharge and controllable load dispatch within the next 1-4 hours based on the rolling updated prediction information, and its optimization objective is to minimize the deviation cost:

[0338]

[0339] The constraint conditions include:

[0340] P grid (t) = P load (t) - P PV，pred (t) - P dis (t) - P ch (t);

[0341] P lod (t) = P ctrl (t) + P non-ctrl，pred (t);

[0342] 0 ≤ P ch (t) ≤ P ch，max ;

[0343] 0 ≤ P dis (t) ≤ P dis，max ;

[0344] P ch (t) · P dis (t) = 0;

[0345]

[0346] SOC min ≤ SOC(t) ≤ SOC max ;

[0347] 0 ≤ P ctrl (t) ≤ P ctrl，max ;

[0348]

[0349] Tmin ≤T in (t) ≤ T max ;

[0350] Among them, C dev (t) represents the penalty cost of grid power deviation (yuan / kWh), P grid，plan (t) represents the grid interaction power of the day-ahead plan (kW), C ctrl (t) represents the penalty cost of controllable load adjustment (yuan / kWh), P ctrl，ref (t) represents the reference power of the controllable load (kW), H represents the optimization time domain length (hours), and Δt represents the time step (hours).

[0351] Considering the comfort constraint of the controllable load, the PMV index is introduced:

[0352] PMV(t) = f(T in (t), RH(t), v air (t), MET, CLO);

[0353] -0.5 ≤ PMV(t) ≤ 0.5;

[0354] Considering the continuous operation time constraint of the controllable load:

[0355] T on (t) ≥ T min，on ·[z(t) - z(t - 1)];

[0356] T off (t) ≥ T min，off ·[z(t - 1) - z(t)];

[0357] Among them, T on (t) and T off (t) respectively represent the duration of turning on and turning off, T min，on and T min，off respectively represent the minimum turning-on and turning-off times, and z(t) represents the device state (1 represents turning on, 0 represents turning off).

[0358] 5.3 Minute-level real-time control layer:

[0359] The minute-level real-time control layer optimizes the charging and discharging power of energy storage and the power of controllable loads based on real-time measurement data and short-term prediction results to achieve real-time balance of energy supply and demand. Its MPC optimization problem is expressed as:

[0360]

[0361] The constraint conditions include:

[0362] P grid(t + i|t) = P load (t + i|t) - P PV (t + i|t) - P dis (t + i|t) + P ch (t + i|t);

[0363] P load (t + i|t) = P ctrl (t + i|t) + P non-ctrl (t + i|t);

[0364] 0 ≤ P ch (t + i|t) ≤ P ch，max ;

[0365] 0 ≤ P dis (t + i|t) ≤ P dis，max ;

[0366] P ch (t + i|t)·P dis (t + i|t) = 0;

[0367]

[0368] SOC min ≤ SOC(t + i|t) ≤ SOC max ;

[0369] 0 ≤ P ctrl (t + i|t) ≤ P ctrl，max ;

[0370] ΔP ctrl，min ≤ ΔP ctrl (t + i|t) ≤ ΔP ctrl，max ;

[0371]

[0372] T min ≤ T in (t + i|t) ≤ T max ;

[0373] Among them, P grid (t + i|t) represents the predicted value of the grid power at time t + i (kW), P grid，ref (t + i) represents the reference grid power (kW), SOC(t + i|t) represents the predicted value of the state of charge at time t + i, SOC ref (t + i) represents the reference state of charge, ΔP ctrl(t + i|t) represents the change in controllable load power (kW), w1, w2, and w3 are weight coefficients, and N represents the prediction time domain length (usually 10 - 15 minutes).

[0374] To improve the computational efficiency, the above - mentioned problem is formulated in the state - space form:

[0375] x(t + 1) = Ax(t)+Bu(t)+Ed(t);

[0376] y(t) = Cx(t);

[0377] where x(t)=[SOC(t), T in (t)] T represents the state vector, u(t)=[P ch (t), P dis (t), P ctrl (t)] T represents the control - input vector, d(t)=[P PV (t), P non-ctrl (t), T out (t)] T represents the disturbance vector, y(t)=[P grid (T), SOC(t), T in (t)] T represents the output vector, and A, B, E, and C are system matrices.

[0378] Considering the charge - discharge mutual - exclusion constraint (P ch (t)·P dis (t)=0), the original problem is transformed into a mixed - integer quadratic programming (MIQP) problem:

[0379] P ch (t)≤P ch，max ·δ(t);

[0380] P dis (t)≤P dis，max ·(1 - δ(t));

[0381] where δ(t)∈{0, 1} is a binary variable, δ(t)=1 represents charging, and δ(t)=0 represents discharging.

[0382] 5.4 Hierarchical Coordination Mechanism:

[0383] To ensure the coordination and consistency between different control levels, the present invention designs the following coordination mechanism:

[0384] 1. Target Constraint Transmission: The optimization results of the upper-layer control serve as the constraints or reference trajectories for the lower-layer control. Specifically, the energy storage SOC trajectory of the day-ahead planning layer serves as the reference trajectory for the hourly scheduling layer; the grid power and SOC trajectory of the hourly scheduling layer serve as the reference trajectories for the minute-level control layer.

[0385] 2. Prediction Information Sharing: The same prediction results are used at different levels, but different time scales and levels of detail are considered. The day-ahead planning layer uses low-frequency trend term prediction; the hourly scheduling layer comprehensively considers low-frequency trend terms and medium-frequency cycle terms for prediction; the minute-level control layer simultaneously considers low-frequency trend terms, medium-frequency cycle terms, and high-frequency fluctuation terms for prediction. The prediction results are stored in a unified data structure form and aggregated or refined according to different time scales.

[0386] 3. Result Feedback Mechanism: The execution results of the lower-layer control are fed back to the upper-layer control as a reference for the upper-layer control to re-plan. For example, the actual energy storage SOC state executed by the minute-level control is regularly fed back to the hourly scheduling layer to correct the hourly scheduling plan; the execution results of the hourly scheduling are fed back to the day-ahead planning layer every few hours to adjust the plan for subsequent periods.

[0387] 4. Time-domain Rolling Optimization: Each layer of control adopts the rolling time-domain optimization method. At the end of each control step, based on the latest system state and prediction information, the optimization problem is re-solved. The expression for rolling optimization is:

[0388]

[0389] U * (t) = [u * (t|t), u * (t + 1|t),..., u * (t + N - 1|t)];

[0390] u applied (t) = u * (t|t);

[0391] where, J * (t) represents the optimal objective function value at time t, U(t) represents the control sequence, L(·) represents the stage cost function, U * (t) represents the optimal control sequence, and u appliied (t) represents the control input actually applied.

[0392] 5. Inter-layer Relaxation Constraints: To avoid infeasible solutions caused by inter-layer objective conflicts, relaxation variables and relaxation constraints are introduced. For example, when the SOC trajectory of the day-ahead planning is used as a reference for the hourly scheduling, a certain range of deviation is allowed:

[0393] SOC ref(t)-ΔSO Callow ≤SOC(t)≤SOC ref (t)+ΔSOC allow ;

[0394] where SOC ref (t) represents the reference SOC value, and ΔSOC allow represents the allowable deviation range.

[0395] 6. Multi-time scale state estimation: To coordinate the state information of different time scales, a multi-time scale Kalman filter is adopted to achieve accurate estimation of state variables. Its model is:

[0396] Prediction step:

[0397]

[0398] Update step:

[0399]

[0400] P k|k =(I-K k H k )P k|K-1 ;

[0401] where represents the prior state estimate, represents the posterior state estimate, P k|k-1 represents the prior covariance matrix, P k|k represents the posterior covariance matrix, K k represents the Kalman gain, F k represents the state transition matrix, B k represents the control matrix, H k represents the observation matrix, Q k represents the process noise covariance, R k represents the observation noise covariance, z k represents the observation vector.

[0402] 5.5 Multi-objective optimization architecture:

[0403] The real-time optimization of the intelligent building green energy system involves multiple conflicting objectives, including economy (minimizing operating costs), reliability (stability of power balance), environmental friendliness (minimizing carbon emissions), and comfort (maximizing user comfort), etc. The present invention uses the weighted sum method to transform the multi-objective problem into a single-objective problem:

[0404]

[0405] where J represents the comprehensive objective function, J i(x, u) represents the \(i\)-th sub-objective function, and \(w\) i represents the weight coefficient of the \(i\)-th sub-objective, and \(m\) represents the number of sub-objectives.

[0406] Since there may be significant differences in the dimensions and orders of magnitude of different sub-objectives, normalization is required:

[0407]

[0408] Among them, represents the normalized sub-objective function, and represent the minimum and maximum values of the sub-objective function, respectively.

[0409] The definitions of each sub-objective function are as follows:

[0410] 1. Economic objective: Minimize the system operation cost, including the cost of purchasing electricity from the power grid, the depreciation cost of energy storage, and the equipment maintenance cost, etc.

[0411]

[0412] 2. Reliability objective: Minimize the power deviation to ensure the stable operation of the system.

[0413]

[0414] 3. Environmental protection objective: Minimize the carbon emissions.

[0415]

[0416] Among them, \(EF\) grid (t) represents the carbon emission factor of the power grid electricity (\(kg CO_2 / kWh\)).

[0417] 4. Comfort objective: Maximize the user comfort, or minimize the discomfort.

[0418]

[0419] Among them, \(PMV(t)\) represents the predicted mean vote.

[0420] Considering the different focuses of each level of control, the selection of weight coefficients also varies: The day-ahead planning layer pays more attention to the economic and environmental protection objectives; the hourly scheduling layer balances the economic, reliability, and comfort; the minute-level control layer mainly focuses on reliability and comfort.

[0421] Table 3 lists the target weight settings for each level of control:

[0422] Control level Economic weight Reliability weight Environmental protection weight Comfort weight Day - ahead planning layer 0.5 0.1 0.3 0.1 Hourly scheduling layer 0.3 0.4 0.1 0.2 Minute - level control layer 0.1 0.6 0.0 0.3 ; 5.6 Safety constraint consideration:

[0423] To ensure the safe and stable operation of the system, the following safety constraints are introduced in this invention:

[0424] 1. Power ramp rate constraint: Restrict the rate of change of the power of each device to avoid large power fluctuations.

[0425] -r down ≤P(t)-P(t - 1)≤r up ;

[0426] where r down and r up represent the maximum downward and upward rates (kW / min), respectively.

[0427] 2. Energy storage state of charge constraint: Restrict the state of charge of the energy storage system within a safe range to avoid overcharging or over-discharging.

[0428] SOC min ≤SOC(t)≤SOC max ;

[0429] The SOC constraints may vary under different operating modes, such as normal mode, emergency mode, etc.

[0430] 3. Equipment start-stop constraint: Restrict the start-stop frequency of the equipment to avoid reducing the equipment life caused by frequent start-stop.

[0431] T on (t)≥T min,on ·[z(t)-z(t - 1)];

[0432] T off (t)≥T min,of f·[z(t - 1)-z(t)];

[0433] where T on (t) and T off (t) represent the duration of start-up and shutdown respectively, T min,on and T min,off represent the minimum start-up and shutdown times respectively, and z(t) represents the equipment state (1 represents start-up, 0 represents shutdown).

[0434] 4. Power factor constraint: Ensure that the system power factor meets the grid requirements and reduce reactive power loss.

[0435] PF min ≤PF(t)≤1;

[0436] where PF(t) represents the power factor at time t, and PF min represents the minimum power factor requirement (usually 0.9 or 0.95).

[0437] 5. Harmonic Constraint: Control the harmonic content of the system to ensure power quality.

[0438] THD(T) ≤ THD max ;

[0439] where THD(t) represents the total harmonic distortion rate, and THD max represents the maximum allowable value (usually 5%).

[0440] 6. Grid Interaction Constraint: Limit the interaction power with the grid within the allowable range.

[0441] P grid，min ≤ P grid (t) ≤ P grid，max ;

[0442] where P grid,min and P grid,max represent the maximum power selling and maximum power purchasing respectively (kW).

[0443] 7. Temperature Comfort Constraint: Ensure the indoor temperature is within the comfortable range.

[0444] T min ≤ T in (t) ≤ T max ;

[0445] where T in (t) represents the indoor temperature (°C), and T min and T max represent the lowest and highest comfortable temperatures respectively.

[0446] 5.7 Fast Solving Algorithm:

[0447] To meet the computational efficiency requirements of real-time control, the present invention designs corresponding fast solving algorithms for optimization problems at different levels:

[0448] 1. Day-ahead Planning Layer: Since it contains integer variables (such as equipment start-stop status), mixed integer linear programming (MILP) is used for solving. To improve the solving efficiency, the following strategies are adopted:

[0449] a. Problem Simplification: Linearize the non-linear constraints, such as approximating the power square term as a piecewise linear function;

[0450] b. Decomposition and Solving: Decompose the 24-hour optimization problem into several sub-problems, such as four 6-hour problems, and then solve them coordinately;

[0451] c. Heuristic Initial Value: Based on the historical optimization results, provide a good initial solution for the solver;

[0452] d. Solving precision control: Set reasonable termination conditions and stop solving in time after obtaining a sufficiently good solution.

[0453] 2. Hourly scheduling layer: mainly for continuous variable optimization problems, solved by quadratic programming (QP) or linear programming (LP). Optimization strategies include:

[0454] a. Preprocessing: Eliminate redundant constraints and reduce the problem scale;

[0455] b. Warm start: Use the optimization results of the previous time period as the initial solution;

[0456] c. Problem scaling: Appropriately scale variables and constraints to improve numerical stability;

[0457] d. Parallel computing: Use multi-core processors to solve the optimization problems of multiple scenarios in parallel.

[0458] 3. Minute-level control layer: Needs to complete the solution within seconds, using explicit MPC or fast gradient method. Implementation strategies include:

[0459] a. Explicit MPC: Calculate the control law offline and only need to look up the table online;

[0460] b. Fast gradient method: Utilize the structural characteristics of the problem to design an efficient iterative algorithm;

[0461] c. Problem simplification: Use a linearized model to reduce the computational complexity;

[0462] d. Computational code optimization: Write the core algorithm in C / C++ and improve the execution efficiency through compilation optimization.

[0463] The basic idea of explicit MPC is to divide the state space into several polyhedral regions and pre-compute the optimal control law for each region:

[0464] u(t) = K i x(t) + g i , if x(t) ∈ R i ;

[0465] where K i and g i represent the control gain matrix and bias vector of the i-th region respectively, and R i represents the i-th state region.

[0466] The iterative formula of the fast gradient method is:

[0467]

[0468] x k+1 = y k+1 + βk (y k+1 -y k );

[0469] where x k and y k represent the variable values at the k-th iteration, Proj Ω represents the projection onto the feasible set Ω, α represents the step size, represents the gradient of the objective function, and β k represents the momentum parameter.

[0470] 6. Adaptive parameter adjustment:

[0471] To improve the adaptability and robustness of the system, the present invention adopts an adaptive parameter adjustment mechanism to adjust the model parameters and control parameters in real time according to the prediction error and control performance.

[0472] 6.1 Adaptive adjustment of prediction model parameters:

[0473] The adaptive adjustment of the prediction model parameters is based on an online learning algorithm, and the update formula is:

[0474]

[0475] where θ (t) represents the model parameters at time t, η represents the learning rate, represents the gradient of the loss function with respect to the parameters.

[0476] The loss function is defined as the weighted sum of squares of the prediction errors:

[0477]

[0478] where P actual (t - i) represents the actual power at time t - i, and P pred (t - i|θ) represents the predicted power at time t - i based on the parameters θ, α i represents the time decay weight, and M represents the length of the historical time window considered.

[0479] The calculation formula for the time decay weight is:

[0480] α i = λ i-1 ;

[0481] where λ ∈ (0, 1) represents the decay factor, with the most recent data having a larger weight and the weight of earlier data gradually decreasing.

[0482] To prevent overfitting, a regularization term is introduced:

[0483] L regJ(θ) = L(θ) + γ·||θ|| 2 ;

[0484] where γ represents the regularization coefficient and ||θ|| 2 represents the L2 norm of the parameters.

[0485] In practical applications, considering computational efficiency, the Stochastic Gradient Descent (SGD) or Mini-batch Gradient Descent (Mini-batch GD) method is adopted:

[0486]

[0487] where L B represents the loss function calculated on the mini-batch data B.

[0488] The learning rate η is selected using an adaptive method, such as the Adam optimizer:

[0489]

[0490]

[0491] where m t and v t represent the first and second moment estimates respectively, β1 and β2 represent the decay rates of the moment estimates, and represent the moment estimates after bias correction, ∈ is a small constant to prevent division by zero errors.

[0492] 6.2 Adaptive adjustment of MPC control parameters:

[0493] The adaptive adjustment of MPC control parameters mainly targets the weight coefficients w1, w2, and w3, as well as the prediction horizon length N.

[0494] The adaptive adjustment of the weight coefficients is based on performance metric evaluation, and the update formula is:

[0495]

[0496] where represents the j-th weight coefficient at time t, γ j represents the adjustment coefficient, represents the corresponding performance metric.

[0497] The performance metric is defined as:

[0498] [[ID=6 5]]

[0499] where K represents the length of the evaluation time window.

[0500] To avoid the weight coefficient being too large or too small, a constraint is introduced:

[0501]

[0502] The weight coefficient needs to be normalized to ensure that their sum is 1:

[0503]

[0504] The adaptive adjustment of the prediction horizon length N is based on the control stability assessment, and the update formula is:

[0505] N (t-1) = max{N min , min{N max , N (t) + δ·sign(S (t) - S threshold )}};

[0506] Among them, N min and N max represent the minimum and maximum allowable prediction horizon lengths respectively, δ represents the adjustment step size, S (t) represents the control stability index, and S threshold represents the stability threshold.

[0507] The control stability index is defined as:

[0508]

[0509] The selection of an appropriate value of N is crucial for the MPC control performance: too small an N will lead to myopic decision-making and poor control performance; too large an N will increase the computational complexity and may introduce more prediction errors. The present invention achieves a balance between control performance and computational efficiency through adaptive adjustment.

[0510] In addition to the weight coefficient and the prediction horizon length, other MPC parameters such as the control horizon length and the constraint softening coefficient can also be adaptively adjusted to make the control system better adapt to environmental changes and uncertainties.

[0511] 6.3 Adaptive adjustment case:

[0512] Taking the photovoltaic power mutation scenario as an example, the action mechanism of the adaptive parameter adjustment is illustrated:

[0513] 1. Under normal conditions, the system operates stably, w1 = 0.6 (power deviation weight), w2 = 0.3 (SOC weight), w3 = 0.1 (control stability weight), N = 10 (prediction horizon length).

[0514] 2. When a significant fluctuation in photovoltaic power is detected (judged by the power change rate or spectral characteristics), the system automatically adjusts the parameters: increase w1 to 0.8, decrease w2 to 0.15, decrease w3 to 0.05, shorten N to 8, and improve the response speed of the system to disturbances.

[0515] 3. After the photovoltaic power fluctuation ends, the system gradually adjusts the parameters back to the normal values according to the actual control performance to ensure long-term stable operation.

[0516] Table 4 shows the adaptive parameter settings under different scenarios:

[0517] Operating scenario Power deviation weight w_1 SOC weight w_2 Control smoothness weight w_3 Prediction horizon length N Normal operation 06 0.3 01 10 Photovoltaic power mutation 0.8 0.15 0.05 8 Load demand mutation 0.7 0.2 01 9 Peak electricity price period 0.5 0.4 0.1 12 Grid fault 0.3 0.6 0.1 6 ;

[0518] The present invention has been experimentally verified in an intelligent office building for 3 months. The building is equipped with a 100 kWp rooftop photovoltaic system and a 200 kWh lithium battery energy storage system. The main loads include lighting, air conditioning, office equipment, etc., and the peak load is about 120 kW.

[0519] The experiment compared three control methods:

[0520] Method A: Traditional rule-based control method;

[0521] Method B: Standard MPC method based on ARIMA prediction;

[0522] Method C: The method of the present invention.

[0523] Experimental design:

[0524] 1. Data collection stage: Collect historical data such as photovoltaic output, load demand, and environmental parameters of the building for model training and parameter adjustment.

[0525] 2. Model training stage: Based on the collected historical data, train the prediction model and the optimal control model to determine the initial parameter settings.

[0526] 3. Control experiment stage: Divide the building into three functionally similar areas and apply the three control methods respectively. Rotate the control areas once a week to eliminate the influence of regional differences.

[0527] 4. Special scenario test: During the experiment, specifically set some special scenarios, such as sudden changes in photovoltaic power (simulated by partial shading), sudden changes in load (by controlling the start and stop of equipment), etc., to test the anti-disturbance ability of the system.

[0528] Evaluation indicators:

[0529] 1. Prediction accuracy: Evaluate the prediction accuracy using indicators such as the mean absolute percentage error (MAPE) and the root mean square error (RMSE).

[0530] 2. Economy: Calculate the daily average operating cost, including the cost of purchasing electricity from the power grid, the depreciation cost of energy storage, etc.

[0531] 3. Reliability: Evaluate the system stability using indicators such as the standard deviation of power fluctuation and the system recovery time.

[0532] 4. Environmental friendliness: Calculate the carbon emissions based on the carbon emission factor of the grid electricity and the actual electricity consumption.

[0533] 5. Comfort: Evaluate the user comfort using the PMV index and the temperature deviation.

[0534] 6. Comprehensive performance: Comprehensively consider the above indicators and calculate the weighted performance score.

[0535] The experimental results are shown in Table 5:

[0536] Performance index Method A Method B Method C (This invention) Photovoltaic prediction mean absolute percentage error (MAPE) 15.2％ 8.7％ 6.3％ Load prediction mean absolute percentage error (MAPE) 12.4％ 7.5％ 5.8％ Photovoltaic local consumption rate 68.5％ 76.3％ 83.7％ Daily average operating cost (yuan) 426.8 389.5 347.2 Standard deviation of power fluctuation (kW) 12.3 8.6 5.4 System recovery time (minutes) 8.5 5.2 3.1 Carbon emission (kg CO / day) 278.6 245.3 213.6 Temperature comfort deviation (℃) 1.8 1.2 0.7 Mean absolute value of PMV 0.75 0.53 0.31 Comprehensive performance score (out of 100) 62.5 78.3 91.7 ;

[0537] The data in Table 5 show that compared with the traditional method, the present invention has significant improvements in various performance indicators:

[0538] The present invention provides a theoretical support and a technical solution for the efficient operation of the green energy system of intelligent buildings, effectively solving the problems of strong time-variability, low prediction accuracy, and lagging control response existing in the prior art, and has a certain significance for promoting building energy conservation and emission reduction and achieving the "dual carbon" goal.

Claims

1. An optimized power management method for an intelligent building green energy system, characterized in that, Including the following steps: Establish a mathematical model of an intelligent building energy system including a photovoltaic system, an energy storage system, a controllable load, and grid interaction; Use short-time Fourier transform to perform time-frequency domain analysis on photovoltaic output and load demand, extract frequency domain characteristics of multiple time scales, and decompose the time series into a trend term, a periodic term, and a random fluctuation term; Based on the Fourier time-frequency characteristics, construct a multi-time scale prediction model for photovoltaic output and load demand; Establish a hierarchical MPC control framework including a day-ahead planning layer, an hourly economic dispatch layer, and a minute-level real-time control layer; Achieve real-time balance of energy supply and demand at the minute-level real-time control layer; Adjust the model parameters and control parameters through an online learning algorithm.

2. The method according to claim 1, characterized in that, The calculation formula for the output power of the photovoltaic system is: P PV (t) = η PV ·S·G(t)·[1 - β(T c (t) - T ref )]; Among them, P PV (t) represents the photovoltaic output power at time t, η PV represents the conversion efficiency of the photovoltaic module, S represents the total area of the photovoltaic module, G(t) represents the solar radiation intensity at time t, β represents the temperature coefficient, T c (t) represents the module temperature at time t, T ref represents the reference temperature under standard test conditions.

3. The method according to claim 1, characterized in that The dynamic model of the energy storage system is expressed as: Among them, SOC(t) represents the state of charge of the energy storage at time t, η ch represents the charging efficiency, η dis represents the discharging efficiency, P ch (t) represents the charging power at time t, P dis (t) represents the discharging power at time t, Δt represents the time step, E cap represents the rated capacity of the energy storage system.

4. The method according to claim 1, wherein The controllable load is a thermoelectric device, and its model is: Among them, T in (t) represents the indoor temperature at time t, C th represents the heat capacity of the building, Q env (t) represents the environmental heat exchange power at time t, and COP represents the energy efficiency coefficient of the device; Q env (t) = UA·(T out (t) - T in (t)) + Q int (t) + Q sol (t); Among them, UA represents the total heat transfer coefficient of the building, and T out (t) represents the outdoor temperature, and Q int (t) represents the internal heat gain, and Q sol (t) represents the solar radiation heat gain.

Images