A hierarchical distributed control method and system for power grid resilience

Abstract

This invention proposes a hierarchical distributed control method and system for power grid resilience, belonging to the field of power system operation and control. The method includes: estimating the frequency-voltage joint sensitivity of active and reactive power disturbances at each node based on real-time power system operation data and power network topology, obtaining a sensitivity tensor and dynamically partitioning it; at the partition granularity, combining the timing constraints of relay protection, low-frequency and low-voltage load shedding, reactive power compensation switching, and distribution switch operations, outputting a coordinated action sequence; publishing scenario parameters and consistency conditions through the master station, with each partition's edge controller solving local subproblems in parallel and performing consistency verification; implementing power-energy dual-dimensional allocation and reserve capacity reservation for energy storage within the partition, with upper-layer event-driven setting of quotas and safety domains for energy storage state of charge, and lower-layer providing virtual inertia support, primary frequency regulation, and ramp compensation. This invention effectively improves the disturbance resilience of novel power systems.

Description

Technical Field

[0001] This invention belongs to the field of power system operation and control technology, and in particular relates to a hierarchical distributed control method and system for power grid resilience. Background Technology

[0002] The statements in this section are merely background information related to the present invention and do not necessarily constitute prior art.

[0003] With the large-scale grid connection of renewable energy and the widespread penetration of power electronic devices, modern power systems exhibit typical characteristics of low inertia, weak damping, and strong coupling, with significantly enhanced spatiotemporal coupling characteristics of frequency and voltage. Against this backdrop, power systems face more severe disturbance challenges; a single fault can trigger a chain reaction, leading to widespread power outages. Enhancing the resilience of the power grid under extreme disturbance scenarios has become a core concern in the field of power system operation and control.

[0004] Currently, emergency control of power systems mainly relies on traditional methods such as low-frequency load shedding (LFLS), low-voltage load shedding (UVLS), relay protection, and reactive power compensation. However, existing technical solutions have the following prominent problems: First, both LFLS and UVLS generally employ fixed action thresholds and static partitioning strategies, lacking dynamic adaptability to changes in system operating status and the spatiotemporal distribution of adjustable resources. When significant changes occur in system inertia levels, load composition, or distributed power output, fixed settings are difficult to match actual operating requirements, easily leading to under-shearing or over-shearing and worsening disturbance recovery performance.

[0005] Second, relay protection, low-frequency and low-voltage load shedding, reactive power compensation switching, and distribution switch operation operate at different time scales, and existing methods lack a unified time-series coordination framework. Each control link is set independently and operates separately, which can easily lead to problems such as timing misalignment, action conflict, or repeated disconnection, resulting in the overall effectiveness of the control strategy being compromised.

[0006] Third, in the context of highly permeable distributed adjustable resources, centralized optimization methods rely on global information collection and central decision-making, resulting in heavy communication burdens and long computation times, making it difficult to meet the timeliness and reliability requirements of control in disturbed scenarios. Furthermore, centralized solutions are not conducive to protecting the privacy of operations in different areas, face the risk of single points of failure in the event of communication anomalies, and lack system robustness. Summary of the Invention

[0007] To overcome the shortcomings of the prior art, this invention provides a hierarchical distributed control method and system for grid resilience, aiming to solve the technical problems of static partitioning lacking adaptability, lack of unified coordination across time scales among multiple types of control resources, and the temporary effectiveness and reliability bottlenecks of centralized optimization in communication and computing in existing power system emergency control schemes.

[0008] To achieve the above objectives, one or more embodiments of the present invention provide the following technical solutions: The first aspect of this invention provides a hierarchical distributed control method for grid resilience; A hierarchical distributed control method for grid resilience includes: Based on real-time data of power system operation and power network topology, the frequency-voltage joint sensitivity of active power disturbance and reactive power disturbance is estimated at each node, and the sensitivity tensor is obtained and dynamically partitioned. By combining the timing constraints of relay protection, low-frequency and low-voltage load shedding, reactive power compensation switching and distribution switch operation at the partition granularity, a cross-time scale collaborative optimization model is established to output a coordinated action sequence. The scene parameters and consistency conditions are published by the main station, and the edge controllers of each partition solve the local sub-problems in parallel, exchanging only the boundary quantities and performing consistency checks. The system implements power-energy dual-dimensional allocation and reserve capacity reservation for energy storage within the zone. The upper layer uses event-driven settings to define the quota and the safety domain of the energy storage charge state, while the lower layer provides virtual inertia support, primary frequency regulation, and ramp compensation.

[0009] As a further technical solution, based on real-time data of power system operation and power network topology, the frequency-voltage joint sensitivity of active power disturbances and reactive power disturbances is estimated at each node, resulting in a sensitivity tensor which is then dynamically partitioned, including: When any phasor measurement unit detects that the frequency change exceeds a preset frequency change threshold, the voltage change exceeds a preset voltage change threshold, or the frequency change rate exceeds a preset frequency change rate threshold, a sliding window of a preset length is activated. Within the sliding window, the sensitivity of each node is estimated using the recursive least squares method, and the sensitivity tensor representation of the node is obtained. The electrical coupling weights between nodes are calculated based on the sensitivity tensor and an electrical coupling graph is constructed. A graph clustering mixed integer programming algorithm with resource coverage constraints is used to complete the dynamic partitioning of the power grid.

[0010] As a further technical solution, the method of calculating the electrical coupling weights between nodes based on the sensitivity tensor and constructing an electrical coupling graph, and then using a graph clustering mixed-integer programming algorithm with resource coverage constraints to complete the dynamic partitioning of the power grid, includes: Using nodes as vertices and lines as edges, an electrical coupling graph is constructed by weighting the edges with the weighted sum of sensitivity spectrum differences, electrical distances, and frequency correlation coefficients: The electrical coupling graph is partitioned using a spectral clustering algorithm with resource coverage constraints. The goal is to maximize the electrical coupling strength within each partition, with constraints that each partition contains at least one fast-response resource and one slow-response resource, and that the power of cross-regional tie lines does not exceed the upper limit. This results in a dynamic partitioning scheme and the boundaries of each partition.

[0011] As a further technical solution, in the process of establishing a cross-timescale collaborative optimization model, a unified discrete time axis is adopted, a basic sampling step size is set, and the minimum hold time of the control action and the minimum time interval between adjacent actions are used to construct the step constraint in order to express the inherent timescale characteristics of different control actions. The objective function of the collaborative optimization model is a multi-objective weighted function: in Indicates load level The weighting coefficients, Indicates the equivalent load shedding power. This indicates a penalty for exceeding the frequency deviation limit. Indicates the frequency deviation safety threshold. This indicates the penalty for voltage deviation exceeding the limit. Represents a node Voltage deviation safety threshold, Indicates switch operation indication. Indicates the cost of switching operations. , and These are the penalty coefficients for the corresponding items.

[0012] As a further technical solution, scene parameters and consistency conditions are published through the main station, and each partition edge controller solves local sub-problems in parallel, exchanging only boundary quantities and performing consistency checks, including: The main station uniformly distributes the disturbance scenario set, global security threshold, and boundary consistency constraints to all partition edge controllers. Each zone independently constructs a local sub-optimization model with energy storage, reactive power, controllable load, and switch as optimization variables, and only injects power, node phase angle, and node voltage boundary quantities into the boundary of the interconnection line between adjacent zones. The alternating direction multiplier method is used to complete local optimization, boundary consensus averaging, and dual multiplier iterative update. The global boundary consistency check is completed based on the dual convergence criteria of the original residual and the dual residual.

[0013] As a further technical solution, the upper-layer event-driven setting of quotas and energy storage state of charge security domains includes: Determined in the event Period, Zone The maximum available power that energy storage can provide instantaneously. and the cumulative available energy For each event With partitions The energy control law for BESS energy storage is as follows: in, express The energy of the energy storage device at all times Indicates the energy storage charging power. Indicates the energy storage discharge power, parameters and These represent the charging and discharging efficiencies of energy storage, respectively. Indicates the energy storage reserve capacity. Indicates energy storage and backup energy. Indicates the minimum state of charge. Indicates the maximum state of charge. This indicates the rated storage capacity of the energy storage system.

[0014] As a further technical solution, the lower layer provides virtual inertia support, primary frequency modulation, and ramp compensation, including: The lower-level fast loop output energy storage power control law is as follows, and ramp compensation is performed within a certain time to ensure that the net imbalance power of the zone is continuously transferred to the upper-level plan: in, For partitioning Energy storage Active power command at any given time; The virtual inertia coefficient for energy storage; For partitioning The rate of change of frequency; This refers to the primary frequency regulation ratio for energy storage. For partitioning Frequency deviation; The integral coefficient for energy storage frequency recovery; This represents the maximum available active power capacity for energy storage. This is the sampling step size.

[0015] A second aspect of the present invention provides a hierarchical distributed control system for grid resilience.

[0016] A hierarchical distributed control system for grid resilience includes: The main site is used to publish scene parameters, global control thresholds, and consistency conditions required for alternating direction multiplier method iteration; Multiple partition edge controllers are configured, with one partition edge controller for each partition. Each partition edge controller has a built-in deterministic real-time kernel, a state estimation module, a local optimization module, and a boundary coordination module. The state estimation module is used to estimate the system state within the partition based on measurement data. The local optimization module is used to solve the local sub-problems corresponding to the partition. The boundary coordination module is used to exchange boundary quantity information with the partition edge controllers of adjacent partitions and perform consistency checks. The measurement unit includes phasor measurement units distributed at various nodes of the power system, used to provide synchronous phasor measurement and time synchronization signals; The execution unit includes an energy storage device, a reactive power compensation device, a controllable load, and a power distribution switchgear, and is used to receive and execute the action commands output by the partition edge controller.

[0017] A third aspect of the present invention provides a computer-readable storage medium having a program stored thereon, which, when executed by a processor, implements the steps of a hierarchical distributed control method for grid resilience as described in the first aspect of the present invention.

[0018] A fourth aspect of the present invention provides an electronic device including a memory, a processor, and a program stored in the memory and executable on the processor, wherein the processor executes the program to implement the steps of a hierarchical distributed control method for grid resilience as described in the first aspect of the present invention.

[0019] The above one or more technical solutions have the following beneficial effects: (1) This invention collects data from each node in real time through a phasor measurement unit and uses a recursive least squares method with a forgetting factor to estimate the frequency-voltage joint sensitivity tensor online, which can accurately characterize the dynamic response characteristics of each node to active / reactive disturbances. An electrical coupling diagram is constructed based on sensitivity spectrum differences, electrical distance, and frequency correlation coefficients, and dynamic partitioning is performed using spectral clustering with resource coverage constraints. This enables the partitioning results to adaptively track changes in system operating status and adjustable resource distribution, overcoming the shortcomings of traditional fixed partitioning strategies in terms of insufficient adaptability when system inertia level, load composition, or distributed power output changes. At the same time, by forcing each partition to contain at least one fast-response resource and one slow-response resource, it ensures that each partition has sufficient adjustment capacity reserves after a disturbance occurs, laying a good physical foundation for subsequent multi-timescale collaborative control.

[0020] (2) This invention integrates relay protection, low-frequency load shedding, low-voltage load shedding, reactive power compensation switching, and distribution switch operation into the same discrete time axis for modeling. By using the minimum holding time and the step constraint of the minimum time interval between adjacent actions to express the inherent time scale characteristics of different control actions, a unified optimization framework spanning milliseconds, seconds, and minutes is constructed. By introducing a mixed integer constraint set such as the "compensation before load shedding" timing priority rule, protection-load shedding mutual exclusion constraint, and the limit on the number of switching operations, the timing misalignment and action conflict between relay protection malfunction, frequent reactive power compensation switching, and load shedding actions are avoided from a mechanism perspective. This effectively reduces the risk of over-switching and incorrect switching caused by the independent setting of each control link in traditional methods, and significantly improves the overall coordination and effectiveness of the control strategy.

[0021] (3) This invention adopts a hierarchical distributed architecture of main station publishing - independent partition solution - boundary information exchange - consistency verification. Each partition only exchanges the injected power, voltage amplitude and phase angle estimates of the boundary nodes, without exposing the internal operating parameters and control strategies of the partition, effectively protecting the operating privacy of each partition. The distributed decomposition of the global problem is achieved by the alternating direction multiplier method, and each partition solves the local subproblems in parallel, which greatly reduces the communication bandwidth requirements and central computing pressure, and significantly shortens the calculation time of the control strategy. At the same time, this invention designs a complete anomaly handling mechanism, which automatically switches to conservative control mode or closed solution mode within the partition when there is measurement anomaly, communication interruption or iteration divergence, avoiding the systemic risk caused by single point failure of the main station in the centralized scheme, and significantly improving the robustness and engineering practicality of the control system.

[0022] (4) This invention breaks through the limitations of traditional fixed charging and discharging strategies for energy storage and proposes a dual-dimensional allocation of power and energy and a two-level control architecture. The upper layer allocates power and energy quotas to the energy storage device in an event-driven manner. The power quota ensures that the energy storage has sufficient instantaneous amplitude support capability at the lowest frequency trough or the deepest voltage drop. The energy quota ensures that the regulation capability will not prematurely exit due to insufficient state of charge during the duration of the event, thus achieving a synergistic balance between instantaneous power support and continuous energy guarantee for energy storage. The lower-level fast loop provides rapid power injection to suppress the rate of frequency change at the sub-second level through three-stage coordinated control of virtual inertia response, primary frequency droop and integral compensation. It provides continuous support proportional to the frequency deviation at the second level and completes the smooth power transition within 1 to 5 seconds. This fully taps the rapid regulation potential of energy storage, effectively improves the frequency trough and voltage minimum point of the power grid under severe disturbance scenarios, and reduces the equivalent load shedding.

[0023] Advantages of additional aspects of the invention will be set forth in part in the description which follows, and in part will be obvious from the description, or may be learned by practice of the invention. Attached Figure Description

[0024] The accompanying drawings, which form part of this invention, are used to provide a further understanding of the invention. The illustrative embodiments of the invention and their descriptions are used to explain the invention and do not constitute an improper limitation of the invention.

[0025] Figure 1 This is a flowchart of the method in the first embodiment.

[0026] Figure 2 This is a system structure diagram of the second embodiment. Detailed Implementation

[0027] It should be noted that the following detailed descriptions are exemplary and intended to provide further illustration of the invention. Unless otherwise specified, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains.

[0028] It should be noted that the terminology used herein is for the purpose of describing particular implementations only and is not intended to limit the exemplary implementations of the present invention.

[0029] Where there is no conflict, the embodiments and features in the embodiments of the present invention can be combined with each other.

[0030] To address the challenges of static partitioning lacking adaptability, insufficient coordination of multiple control resources across time scales, and low reliability of centralized optimization in existing power system emergency control systems, this invention first performs dynamic partitioning with adjustable resource constraints based on online estimation of frequency-voltage joint sensitivity. Then, it establishes a unified time-series collaborative optimization model spanning milliseconds, seconds, and minutes at the partition granularity and outputs action sequences. Next, it employs the alternating direction multiplier method to achieve distributed optimization involving master station deployment, independent partition solution, and boundary consistency verification. Finally, it implements power-energy dual-dimensional allocation and two-level control with energy storage as the core. Through these technical approaches, this invention effectively reduces equivalent load shedding, suppresses cross-regional power fluctuations, and significantly improves grid resilience while ensuring frequency / voltage safety.

[0031] Example 1 This embodiment discloses a hierarchical distributed control method for grid resilience. It achieves dynamic partitioning by constructing a frequency-voltage joint sensitivity tensor, establishes a collaborative optimization model across millisecond-second-minute time scales, and uses the alternating direction multiplier method to achieve hierarchical distributed solution. While satisfying frequency and voltage safety constraints, it reduces the equivalent load shedding and recovery time, and suppresses cross-regional power fluctuations.

[0032] like Figure 1 As shown, a hierarchical distributed control method for grid resilience includes: Step S1: Based on real-time data of power system operation and power network topology, estimate the frequency-voltage joint sensitivity of active power disturbance and reactive power disturbance by node, obtain the sensitivity tensor and perform dynamic partitioning.

[0033] Step S11: When any synchronous phasor measurement unit (PMU) detects a frequency deviation exceeding a preset threshold... Voltage deviation exceeds preset threshold Or, if ROCOF exceeds the limit, the startup length is... A sliding window. Among them, the parameters... Indicates the frequency change threshold. This represents the voltage change threshold. The sliding window duration can be adaptively adjusted according to the grid scale; for example, a 200ms sliding window is set for regional distribution networks, and a 500ms sliding window is set for provincial transmission networks. Within the window, frequency, voltage, and active and reactive power time-series data before and after the disturbance are fully collected, filtering out random measurement noise to ensure the accuracy of sensitivity tensor identification. The sliding window adopts a rolling update mechanism, refreshing a set of measurement data every 10ms. A recursive least squares method with a forgetting factor gradually weakens the weight of historical old data, focusing on retaining the dynamic response characteristics at the moment of the disturbance. When the disturbance disappears and the frequency and voltage return to the safe range, the sliding window automatically closes, stopping online sensitivity identification, reducing the computational overhead of the edge controller, and realizing an adaptive computing power scheduling mechanism of "low load under normal conditions and high computing power during disturbances."

[0034] Step S12, for each node and The sliding time after the disturbance is triggered Within the window, the sensitivity vector is estimated using the recursive least squares method. The regression form is: in, This indicates the sensitivity of node active power variations to system frequency. This indicates the coupling effect of node reactive power variations on system frequency. This indicates the local sensitivity of node voltage to active power injection. This indicates the local sensitivity of node voltage to reactive power injection. , , and These are the regression coefficients; The deviation term in the regression formula for frequency change is represented by... This represents the deviation term in the regression formula for node voltage change. , , Representing the injection nodes The active power, reactive power, and node voltage.

[0035] The recursive formula for the least squares method is: in, , This represents the forgetting factor.

[0036] Step S13, with node For vertices and lines For each edge, the weights are the weights of the sensitivity spectrum difference, electrical distance, and frequency correlation coefficient. Construct an electrical coupling diagram The weights of the lines in the network are: in, Represents the set of network lines and , Indicates differences in sensitivity spectrum. Represents a positive definite weight matrix. Indicates electrical distance (the larger the value, the weaker the interaction between voltage and power). This represents the frequency correlation coefficient of the PMU. , , This represents the normalized weight.

[0037] Step S14, in each partition A mixed-integer timing optimization model is established. To ensure that nodes with similar sensitivity, strong electrical coupling, and high dynamic correlation are located in the same partition, and that each partition covers at least one fast resource and one slow resource to reserve some adjustment capability for subsequent control, the following dynamic partitioning MILP model is established: in, Represents a set of partitions, and Discrete sets in the time domain: Sampling step size Minimum resource coverage constraint: in, This indicates a rapid collection of resources, including energy storage, static var compensators (SVC / STATCOM), etc. , This indicates a set of slow resources, including primary frequency regulation of conventional units, interruptible loads, etc. .

[0038] Inter-regional connection line restrictions: The above constraints can be added to the objective function: The adjustable resource set includes the energy storage resource set. Reactive power compensation resource set interruptible load resource set Conventional unit resource set Each partition Must contain at least one With one and satisfy the upper limit of boundary exchange. .

[0039] Spectral clustering combined with small-scale MILP can be used to perform coverage correction to obtain partitions. With boundary .parameter For binary variables, if the node Belongs to partition but ,parameter For auxiliary binary variables, The upper limit of boundary exchanges allowed by the link line. This indicates that it is allowed to put the edge The degree of disconnection, parameters This represents the penalty coefficient for boundary breaks.

[0040] Step S2: At the partition granularity, combine the timing constraints of relay protection, low-frequency low-voltage load shedding, reactive power compensation switching and distribution switch operation to establish a cross-time scale collaborative optimization model and output a coordinated action sequence.

[0041] This step addresses the traditional problem of timing discrepancies among multiple devices by mapping all control actions to a unified discrete time axis, and setting a uniform basic sampling step size to match the PMU synchronous acquisition frequency. The model also incorporates multi-scenario opportunity constraints, considering random fluctuations in renewable energy output, load forecasting errors, and uncertainties in equipment regulation capacity. It employs the Sample Approximation Method (SAA) combined with Conditional Value at Risk (CVaR) to quantify extreme disturbance risks. During optimization, it balances optimal control cost with extreme fault safety margins, avoiding optimization results that only apply to a single standard disturbance scenario and improving the generalization ability of the control strategy.

[0042] To accommodate milliseconds (protection), seconds (compensation, initial load shedding), and minutes (recovery / reconstruction), a unified discrete time axis is adopted. Set the base step size The time scale of different actions is expressed using minimum hold time / minimum interval step constraints. The following objective function (multi-objective weighted) is established: in Indicates load level The weighting coefficients, Indicates the equivalent load shedding power. This indicates a penalty for exceeding the frequency deviation limit. Indicates the frequency deviation safety threshold. This indicates the penalty for voltage deviation exceeding the limit. Represents a node Voltage deviation safety threshold, Indicates switch operation indication. Indicates the cost of switching operations. , and These are the penalty coefficients for the corresponding items.

[0043] A small-signal model of frequency deviation and voltage deviation is constructed. This model represents the factors affecting frequency variation, thus allowing it to be correlated with the objective function to determine which factors influence the frequency deviation and whether it exceeds the limit. As shown below: in, , Representing partitions At any moment Frequency variation and nodes Voltage change Represents the equivalent inertia constant. Represents the load characteristic coefficient. Indicates an active disturbance. Indicates the load shearing amount of LFLS / UVLS. This indicates the injected power for energy storage.

[0044] To ensure the timing coordination and mutual exclusion of control actions across time scales, and to avoid protection malfunctions, overload shedding, and resource conflicts, this embodiment constructs the following timing logic and interlocking constraint system under a unified discrete time axis framework: (1) Trigger threshold and delay constraint (taking LFLS as an example) Let the first The LFLS frequency threshold is The delay threshold is Establish the following constraint functions: in, This represents the perturbation frequency coefficient. Indicates the timer value. Indicates the duration of the action. Indicates the first Duration of level-reduction load, Indicates the first The level of load reduction binary state is engaged or maintained in the load reduction state, and its value is 1.

[0045] (2) Timing interlock constraint for compensation followed by load reduction in, Indicates reactive power compensator or capacitor bank resources The closed state, This indicates that priority will be given to compensating for the delay after the investment. This represents the binary state of the low-frequency load shearing segment switch under the minimum interval time. This indicates the minimum interval time for low-frequency load shedding operations.

[0046] Reactive power compensation hold time constraint: in, This indicates the minimum holding time for reactive power compensation.

[0047] (3) Protection-load reduction mutual exclusion and switching frequency constraint in, Indicates protective element Protective action instructions, When the protection element activates, its state is 1, at which point low-frequency load shearing occurs. Operating status and low-pressure load reduction The action states are all 0. Action indication that shows the binary state of the switching / reconfiguration operation of the protection element.

[0048] Number of switching operations and minimum interval: in, This indicates the maximum number of operations allowed for the protective switch. The minimum duration of a switch operation.

[0049] (4) Load shedding and reactive power equivalent constraints: in, Indicates the low-frequency load shedding power limit. Indicates the low-voltage load shedding power limit (5) Opportunity Constraint (SAA) In the context of the scene Regarding factors such as disturbance magnitude, inertia, and load uncertainty, the probability of frequency / voltage exceeding limits should not exceed 1- This means that the proportion of scenarios that do not violate the constraints must be no less than the specified confidence level. : SAA is converted to .

[0050] It can be linearized using auxiliary binary and the Big M method, or approximated by CVaR. The above formula represents the tail mean of constraint violation severity using CVaR, which limits the worst-case scenario. The severity of violations at the rear has been reduced to a sufficiently small level, among which, Indicates the frequency deviation safety threshold. Indicates the safe threshold for voltage deviation. Represents a scene set, This indicates the confidence level that the scenario opportunity constraint is satisfied. This represents the auxiliary slack variable used in calculating the conditional value at risk (CVaR). This represents the weighted value of the frequency / voltage deviation, i.e., the overall operating cost.

[0051] Step S3: Publish scene parameters and consistency conditions through the main station. Each partition edge controller solves local sub-problems in parallel, only exchanging boundary quantities and performing consistency checks.

[0052] The ADMM distributed iterative process adapts to the real-time computing capabilities of the edge controller, setting up an adaptive penalty coefficient adjustment mechanism. When the boundary residuals are large, the consistency penalty weight is automatically increased to accelerate the convergence speed of the partition boundary; when the residuals gradually decrease, the penalty coefficient is reduced to avoid the optimization target from excessively deviating from the local optimal solution. The global convergence criterion sets dual thresholds: the original residual controls the deviation of power and voltage on both sides of the tie line, and the dual residual controls the deviation of the partition optimization target. The iteration is terminated and control commands are output only when both conditions are met simultaneously. If convergence is not achieved after reaching the maximum number of iterations, the edge controller automatically retrieves the local conservative control plan, relying on the energy storage and controllable loads within the zone to quickly smooth out disturbances, without waiting for the global consistency result, ensuring that control actions are issued in real time and preventing grid instability caused by the lack of command windows.

[0053] Global issues by partition Decomposition, boundary consistency variables are (Tie-up injection, phase angle / voltage) The global problem is decomposed into subproblems in each region. Each region only exchanges the boundary injection power, phase angle, and voltage estimates. The objective function is... It is a partition The operational objectives within this area should be to ensure minimal load shedding, minimal protection actions, and minimal parameter exceedances. The following objective function should be solved: ADMM iteration: <1> Local step: ; <2> Consensus Step: (Calculate the average or weighted average based on the boundary lines); <3> Multiplication step: .

[0054] In the local steps of ADMM, This approach ensures that states deviating from the common boundary are penalized and suppressed, thereby driving the partition boundaries to gradually converge. In the consensus step, the boundaries provided by each partition are weighted and averaged to obtain a new common boundary. In the multiplier step, if there is still a difference between the local boundary and the common state, the duality of the corresponding component is increased, thereby promoting the convergence of communication parameters between adjacent regions. In the local step, each partition independently calculates its optimal control action based on its local state, equipment capabilities, and the previous round's common boundary value. After the local step, it outputs its internal control quantity and uploads the collaborative boundary quantity. In the consensus step, the boundary quantities uploaded by each partition are merged into a new common boundary state. In the multiplier step, the deviation between the local boundary quantity and the common boundary state quantity is recorded, and the deviation is fed back into the next round of local optimization.

[0055] in, Includes , Consistency variables for boundary injection / phase angle / voltage. Indicates the status of the load reduction operation. This indicates the closing status of the reactive power compensator or capacitor bank resources. Indicates the operating status of the protective element. Indicates the energy storage charging and discharging power. It indicates energy storage capacity or energy. This represents the consistency penalty parameter for ADMM adaptive adjustment. Indicates the scaling dual parameter. This represents the boundary selection matrix.

[0056] Convergence criterion: in, Represents the original residual. Represents the dual residual, and represents Original residual convergence threshold This represents the convergence threshold of the dual residual.

[0057] Step S4: Implement power-energy dual-dimensional allocation and reserve capacity reservation for energy storage within the partition. The upper-layer event-driven setting of quotas and safety domains for energy storage charge state provides virtual inertia support, primary frequency regulation, and ramp compensation.

[0058] The upper-layer event-driven module incorporates fault type identification logic. It differentiates between four typical disturbances—short-circuit faults, unit tripping, load surges, and renewable energy disconnection—based on disturbance timing characteristics collected by the PMU. Energy storage power and quotas are allocated differently for each fault type. For example, large-scale renewable energy disconnection is a prolonged active power deficit fault, so the energy storage quota is significantly increased, allowing for sufficient discharge time. Instantaneous line short circuits are short-term impact disturbances, so energy constraints are relaxed to maximize instantaneous power support. Simultaneously, it reserves energy storage standby capacity in real-time, providing adjustment margins to cope with worsening disturbances, and sets upper and lower SOC safety limits to prevent overcharging and over-discharging damage to the energy storage equipment. The lower-layer three-loop control outputs continuous and smooth power commands, eliminating step power surges. A ramp-up transition mechanism ensures seamless integration between the upper-layer optimized long-term power plan and the lower-layer rapid instantaneous support, balancing transient stability and continuous power balance.

[0059] The upper layer implements power and energy control, first determining the event... Period, Zone The maximum available power that energy storage can instantly provide. (Power quota) and cumulative available energy (Energy quota), i.e. For each event With partitions The energy control law for BESS energy storage is as follows: Power quota To ensure peak capacity, energy storage must have sufficient instantaneous amplitude (including virtual inertia / primary frequency regulation / second-level ramp-up) to withstand peak disturbances or reduce ROCOF during periods of lowest frequency or deepest voltage drops; energy quotas. Ensure continuous capability to guarantee that the duration of regulation capability is not too short due to insufficient state of charge SoC during the typical duration of an event.

[0060] The lower-level fast loop output energy storage power control law is as follows, and ramp compensation is performed within 1–5 seconds to ensure that the net imbalance power of the zone is continuously transferred to the upper-level plan: in, express The energy of the energy storage device at all times Indicates the energy storage charging power. Indicates the energy storage discharge power, parameters and These represent the charging and discharging efficiencies of energy storage, respectively. Indicates the energy storage reserve capacity. Indicates energy storage and backup energy. Indicates the minimum state of charge. Indicates the maximum state of charge. This indicates the rated storage capacity of the energy storage system.

[0061] Example 2 This embodiment discloses a hierarchical distributed control system for grid resilience; like Figure 2 As shown, a hierarchical distributed control system for grid resilience includes: A hierarchical distributed control system for grid resilience includes: The main site is used to publish scene parameters, global control thresholds, and consistency conditions required for alternating direction multiplier method iteration.

[0062] As the highest management layer of the system, the master station is responsible for publishing scene parameters, global control thresholds, and consistency conditions required for ADMM iteration. The master station consists of a monitoring host, an engineering workstation, and a communication front-end, employing a dual-machine hot standby configuration to ensure that core equipment has no single point of failure.

[0063] In practice, the master station first generates a shared set of scenario parameters for each sub-region based on the current power grid operating status and disturbance scenario information. These scenario parameters include: disturbance type (such as line short circuit, generator tripping, load surge, etc.), disturbance occurrence time, disturbance duration, and initial system operating status. Global control thresholds include frequency deviation safety thresholds, voltage deviation safety thresholds, and frequency change rate safety thresholds.

[0064] The consistency condition parameters required for ADMM iteration include: the initial value of the consistency penalty parameter, the original residual convergence threshold, the dual residual convergence threshold, and the maximum number of iterations. The master station sends these parameters to each partition edge controller via the communication network and continuously monitors the status information returned by each partition.

[0065] Multiple partition edge controllers are provided, with one partition edge controller for each partition. Each partition edge controller has a built-in deterministic real-time kernel, a state estimation module, a local optimization module, and a boundary coordination module. The state estimation module is used to estimate the system state within the partition based on measurement data. The local optimization module is used to solve the local sub-problems corresponding to the partition. The boundary coordination module is used to exchange boundary quantity information with the partition edge controllers of adjacent partitions and perform consistency checks.

[0066] The deterministic real-time kernel runs a real-time operating system, responsible for task scheduling and timing management, ensuring that each functional module executes under strict time constraints. Specifically, the real-time kernel divides the control cycle into several time slices, including data acquisition time slices, state estimation time slices, local optimization time slices, and boundary coordination time slices. Taking a disturbance-triggered control scenario as an example, the data acquisition time slice is set to 10 milliseconds, the state estimation time slice to 20 milliseconds, the local optimization time slice to 100 milliseconds, and the boundary coordination time slice to 50 milliseconds. Each time slice is executed in a fixed priority order.

[0067] The state estimation module is used to estimate the system state within the partition based on measurement data. Specifically, the state estimation module receives voltage amplitude, voltage phase angle, frequency, active power, and reactive power data from each node of the measurement unit. Weighted least squares method is used for state estimation, with state variables including the voltage amplitude and phase angle of each node. When some measurement data is missing or abnormal, the state estimation module automatically switches to a predictive estimation mode based on historical data and a preset system model, ensuring the continuity and reliability of the state estimation results.

[0068] For the estimation of frequency state variables, the state estimation module uses a recursive least squares method with a forgetting factor to identify the sensitivity tensor of each node online within the sliding time window after the disturbance is triggered, providing the local optimization module with a real-time estimate of the system's dynamic response characteristics.

[0069] The local optimization module is used to solve the local subproblems corresponding to the current partition. In specific implementation, the local optimization module receives the scene parameters and consistency conditions issued by the main station, receives the real-time state estimate of the current partition from the state estimation module, and uses the coordinated action sequence output in step B as the initial iterative value of the local decision variables.

[0070] The objective function of the local subproblem solved by the local optimization module is to minimize the sum of load shedding, frequency / voltage over-limit penalties, and switching operation costs within the local partition. Decision variables include the power shedding at each load level in each time period, the switching status of reactive power compensation devices, the operating status of each segment of low-frequency load shedding, the operating status of each segment of low-voltage load shedding, the operating status of protection components, and the charging and discharging power of energy storage devices.

[0071] The constraints include power balance constraints within this partition, node voltage safety constraints, line power flow constraints, and timing constraints of various control resources. Timing constraints are incorporated into the local optimization model in the form of a mixed-integer linear constraint set, specifically including relay protection action delay and interlocking logic, threshold values ​​and minimum intervals for graded and segmented triggering of low-frequency load shedding, voltage thresholds and delay coordination for low-voltage load shedding, dead zones and switching holding times of reactive power compensation devices, the timing priority rule of "compensation before load shedding," and the maximum number of operations and minimum intervals between adjacent operations for distribution switches.

[0072] The local optimization module uses a mixed-integer linear programming solver. In each step of the ADMM iteration, the local optimization module independently solves the local subproblem of its own partition and passes the solution results (including the injected active power, injected reactive power, voltage magnitude and phase angle estimates at the boundary nodes of the partition) to the boundary coordination module.

[0073] The boundary coordination module is used to exchange boundary quantity information with the partition edge controllers of adjacent partitions and perform consistency checks. In specific implementation, the boundary coordination module establishes a point-to-point communication connection with the boundary coordination modules of adjacent partitions through a communication network. In each ADMM iteration, the boundary coordination module first packages and sends the state information of the boundary nodes of its own partition to the adjacent partitions, while simultaneously receiving the boundary state information sent by the adjacent partitions.

[0074] Upon receiving boundary information from adjacent partitions, the boundary coordination module performs a weighted average of the state values ​​on both sides of the same boundary to obtain a consistent common boundary value, and updates the dual variable of this partition accordingly. The boundary coordination module simultaneously calculates the original residual and the dual residual. If both are less than a preset convergence threshold, it notifies the local optimization module to terminate the iteration and output the final optimization result; otherwise, it continues to the next iteration.

[0075] If the number of iterations exceeds the preset maximum number of iterations and still fails to converge, the boundary coordination module triggers an exception handling mechanism: it issues a "switch to conservative mode" instruction to the local optimization module of this partition, and the local optimization module uses the coordination action sequence output in step B as the final execution instruction.

[0076] The measurement unit, including phasor measurement units distributed at various nodes of the power system, is used to provide synchronous phasor measurements and time synchronization signals.

[0077] The measurement unit includes phasor measurement units (PMUs) distributed at various nodes of the power system, used to provide synchronous phasor measurements and time synchronization signals. Each PMU uses the second pulse signal of the Global Positioning System (GPS) or the BeiDou Navigation Satellite System (BDS) as the synchronization reference to achieve synchronous acquisition of frequency, voltage amplitude, voltage phase angle, and power.

[0078] The measurement unit's data acquisition frequency is no less than 50Hz, with frequency measurement accuracy better than 0.01Hz, voltage phase angle measurement accuracy better than 0.05 degrees, voltage amplitude measurement accuracy better than 0.2%, and power measurement accuracy better than 0.5%. Each PMU uploads the acquired data with precise time stamps to the partition edge controller of its respective partition via the communication network.

[0079] The execution unit includes an energy storage device, a reactive power compensation device, a controllable load, and a power distribution switchgear, and is used to receive and execute the action commands output by the partition edge controller.

[0080] The energy storage device includes a lithium-ion battery energy storage system, a flow battery energy storage system, or a supercapacitor energy storage system. Each energy storage device is equipped with a local controller. The local controller receives power commands issued by the partition edge controller and converts them into pulse width modulation signals of the energy storage converter to realize closed-loop control of the energy storage charging and discharging power.

[0081] The reactive power compensation device includes parallel capacitor banks, static var compensators (SVCs) and / or static synchronous compensators (STATCOMs). Each device is equipped with a switching control terminal, which receives switching commands from the zone edge controller and then performs closing or opening operations.

[0082] Controllable loads include interruptible loads and adjustable loads. Each controllable load is equipped with a load control terminal, which receives load reduction instructions from the zone edge controller and performs load shedding or power reduction operations.

[0083] The power distribution switchgear includes circuit breakers and load switches. Each switchgear is equipped with an intelligent operation terminal, which receives the opening / closing commands issued by the zone edge controller and then executes the corresponding operations.

[0084] After receiving the instruction, all execution units send back an execution confirmation signal and post-execution status feedback information to the partition edge controller, forming a closed-loop control process of "instruction issuance - execution confirmation - status feedback" to ensure that the control instructions are effectively executed.

[0085] Example 3 The purpose of this embodiment is to provide a computer-readable storage medium.

[0086] A computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the steps in a hierarchical distributed control method for grid resilience as described in Example 1.

[0087] The computer-readable storage medium described in this embodiment can be a non-volatile storage medium, including general-purpose storage hardware in the power industrial control field such as solid-state drives (SSDs), industrial-grade USB flash drives, flash memory chips, ROM read-only memory, and distributed storage disk arrays. It is compatible with various devices such as power grid master station servers, local storage units of partition edge controllers, and embedded storage modules of field measurement and control devices. The storage medium contains complete control program code, with each code module corresponding one-to-one with the method steps in Embodiment 1. It is divided into five independent callable modules: sensitivity tensor identification subroutine, dynamic partitioning mixed integer programming solution subroutine, cross-timescale collaborative optimization subroutine, ADMM boundary consistency iteration subroutine, and energy storage power-energy dual-layer control subroutine. Each module has a standardized data interaction interface, which can be loaded and debugged separately, facilitating segmented verification of the algorithm logic by on-site maintenance personnel. The program has multiple sets of typical disturbance scenario test cases pre-built, covering high-frequency faults in new power systems such as new energy grid disconnection, three-phase short circuits on lines, sudden increases in large-capacity loads, and unit tripping. Maintenance personnel can directly call up the scenarios to complete offline simulation verification without rebuilding the simulation examples. Meanwhile, the storage medium's built-in program includes an exception-tolerant subroutine. When packet loss or distortion occurs in the measured data, it automatically initiates sliding window completion and noise filtering logic. If numerical divergence occurs in the iterative calculation, the program automatically switches to a conservative control branch to output basic frequency and voltage regulation commands, preventing program freezes and control function failure. This storage medium can be batch-programmed to multiple edge controllers, enabling standardized deployment of full-domain partitioned control algorithms. After programming, it supports local offline execution of complete control logic. Even if the device is disconnected from the main station's communication network, it can still independently complete a full set of autonomous control processes, including dynamic partitioning, energy storage collaborative regulation, and emergency load reduction, significantly reducing the construction costs and time required for on-site equipment debugging and algorithm upgrades.

[0088] Example 4 The purpose of this embodiment is to provide an electronic device.

[0089] An electronic device includes a memory, a processor, and a program stored in the memory and executable on the processor. When the processor executes the program, it implements the steps in a hierarchical distributed control method for grid resilience as described in Embodiment 1.

[0090] The electronic equipment in this embodiment is divided into two hardware forms: one is a high-performance computing server for the power grid dispatch master station, and the other is an edge real-time control terminal deployed in each distribution zone. Both types of equipment have hardware architectures adapted to the hierarchical distributed control logic of this invention. The processor uses an industrial-grade multi-core real-time processor. The master station equipment is equipped with multiple high-performance CPU parallel computing cores to support synchronous iteration of ADMM across multiple zones of the entire network and pre-simulation of large-scale disturbance scenarios. The zone edge terminals use ARM multi-core processors with hard real-time kernels to ensure strict controllability of millisecond-level PMU data acquisition, online sensitivity identification, and local optimization solution timing. The memory is partitioned into storage areas. High-speed memory is responsible for real-time measurement data and temporary caching of intermediate variables during iteration. Non-volatile memory stores static configuration information such as control programs, power grid topology parameters, energy storage device parameters, and global safety thresholds. The equipment is equipped with a BeiDou / GPS synchronization module to unify the data timescale across the entire network, ensuring the alignment of phasor measurement timing at each node and eliminating time deviation errors in sensitivity tensor identification. The peripheral communication interfaces integrate fiber optic Ethernet, 5G industrial private network, and RS485 power protocol interface, enabling simultaneous connection to all execution units such as PMU measurement devices, energy storage converters, reactive power compensation terminals, and controllable load switches. The equipment features a built-in hardware watchdog and dual-power redundant power supply module, automatically restarting the control program and loading a local conservative control strategy when the processor is overloaded or the program crashes abnormally. The master station-level electronic equipment adopts a dual-machine hot standby architecture; when the main computing unit fails, the backup equipment seamlessly takes over global parameter distribution and global consistency verification tasks. The equipment is equipped with a visual human-machine interface that can display dynamic partition boundaries, energy storage SOC safety domains, ADMM iterative residuals, and frequency and voltage over-limit alarm information in real time, facilitating real-time monitoring of the entire process of grid disturbance control by dispatchers. It also possesses integrated engineering application capabilities including algorithm computation, local autonomy, status visualization, and fault alarm.

[0091] The steps and methods involved in the apparatuses of Embodiments 2, 3, and 4 above correspond to those in Embodiment 1. For specific implementation details, please refer to the relevant description section of Embodiment 1. The term "computer-readable storage medium" should be understood as a single medium or multiple media including one or more instruction sets; it should also be understood as including any medium capable of storing, encoding, or carrying an instruction set for execution by a processor and enabling the processor to perform any of the methods in this invention.

[0092] Those skilled in the art will understand that the modules or steps of the present invention described above can be implemented using general-purpose computer devices. Optionally, they can be implemented using computer-executable program code, thereby allowing them to be stored in a storage device for execution by a computer device, or they can be fabricated as separate integrated circuit modules, or multiple modules or steps can be fabricated as a single integrated circuit module. The present invention is not limited to any particular combination of hardware and software.

[0093] While the specific embodiments of the present invention have been described above in conjunction with the accompanying drawings, this is not intended to limit the scope of protection of the present invention. Those skilled in the art should understand that various modifications or variations that can be made by those skilled in the art without creative effort based on the technical solutions of the present invention are still within the scope of protection of the present invention.

Claims

1. A hierarchical distributed control method for grid resilience, characterized in that, include: Based on real-time data of power system operation and power network topology, the frequency-voltage joint sensitivity of active power disturbance and reactive power disturbance is estimated at each node, and the sensitivity tensor is obtained and dynamically partitioned. By combining the timing constraints of relay protection, low-frequency and low-voltage load shedding, reactive power compensation switching and distribution switch operation at the partition granularity, a cross-time scale collaborative optimization model is established to output a coordinated action sequence. The scene parameters and consistency conditions are published by the main station, and the edge controllers of each partition solve the local sub-problems in parallel, exchanging only the boundary quantities and performing consistency checks. The system implements power-energy dual-dimensional allocation and reserve capacity reservation for energy storage within the zone. The upper layer uses event-driven settings to define the quota and the safety domain of the energy storage charge state, while the lower layer provides virtual inertia support, primary frequency regulation, and ramp compensation.

2. The hierarchical distributed control method for grid resilience as described in claim 1, characterized in that, Based on real-time data of power system operation and power network topology, the frequency-voltage joint sensitivity of active and reactive power disturbances is estimated at the node level. The sensitivity tensor is then dynamically partitioned, including: When any phasor measurement unit detects that the frequency change exceeds a preset frequency change threshold, the voltage change exceeds a preset voltage change threshold, or the frequency change rate exceeds a preset frequency change rate threshold, a sliding window of a preset length is activated. Within the sliding window, the sensitivity of each node is estimated using the recursive least squares method, and the sensitivity tensor representation of the node is obtained. The electrical coupling weights between nodes are calculated based on the sensitivity tensor and an electrical coupling graph is constructed. A graph clustering mixed integer programming algorithm with resource coverage constraints is used to complete the dynamic partitioning of the power grid.

3. The hierarchical distributed control method for grid resilience as described in claim 2, characterized in that, The process involves calculating the electrical coupling weights between nodes based on sensitivity tensors and constructing an electrical coupling graph. A graph clustering mixed-integer programming algorithm with resource coverage constraints is then used to complete the dynamic partitioning of the power grid, including: Using nodes as vertices and lines as edges, an electrical coupling graph is constructed by weighting the edges with the weighted sum of sensitivity spectrum differences, electrical distances, and frequency correlation coefficients: The electrical coupling graph is partitioned using a spectral clustering algorithm with resource coverage constraints. The goal is to maximize the electrical coupling strength within each partition, with constraints that each partition contains at least one fast-response resource and one slow-response resource, and that the power of cross-regional tie lines does not exceed the upper limit. This results in a dynamic partitioning scheme and the boundaries of each partition.

4. The hierarchical distributed control method for grid resilience as described in claim 1, characterized in that, In the process of establishing a collaborative optimization model across time scales, a unified discrete time axis is adopted, a basic sampling step size is set, and the minimum hold time of the control action and the minimum time interval between adjacent actions are used to construct step constraints to express the inherent time scale characteristics of different control actions. The objective function of the collaborative optimization model is a multi-objective weighted function: in Indicates load level The weighting coefficients, Indicates the equivalent load shedding power. This indicates a penalty for exceeding the frequency deviation limit. Indicates the frequency deviation safety threshold. This indicates the penalty for voltage deviation exceeding the limit. Represents a node Voltage deviation safety threshold, Indicates switch operation indication. Indicates the cost of switching operations. , and These are the penalty coefficients for the corresponding items.

5. The hierarchical distributed control method for grid resilience as described in claim 1, characterized in that, The scene parameters and consistency conditions are published through the main site. Each partition edge controller solves local subproblems in parallel, exchanging only boundary quantities and performing consistency checks, including: The main station uniformly distributes the disturbance scenario set, global security threshold, and boundary consistency constraints to all partition edge controllers. Each zone independently constructs a local sub-optimization model with energy storage, reactive power, controllable load, and switch as optimization variables, and only injects power, node phase angle, and node voltage boundary quantities into the boundary of the interconnection line between adjacent zones. The alternating direction multiplier method is used to complete local optimization, boundary consensus averaging, and dual multiplier iterative update. The global boundary consistency check is completed based on the dual convergence criteria of the original residual and the dual residual.

6. The hierarchical distributed control method for grid resilience as described in claim 1, characterized in that, Upper-layer event-driven settings define the security domains for quotas and energy storage state of charge, including: Determined in the event Period, Zone The maximum available power that energy storage can provide instantaneously And the cumulative available energy For each event With partitions The energy control law for BESS energy storage is as follows: in, express The energy of the energy storage device at all times Indicates the energy storage charging power. Indicates the energy storage discharge power, parameters and These represent the charging and discharging efficiencies of energy storage, respectively. Indicates the energy storage reserve capacity. Indicates energy storage and backup energy. Indicates the minimum state of charge. Indicates the maximum state of charge. This indicates the rated storage capacity of the energy storage system.

7. The hierarchical distributed control method for grid resilience as described in claim 1, characterized in that, The lower layer provides virtual inertia support, primary frequency modulation, and ramp compensation, including: The lower-level fast loop output energy storage power control law is as follows, and ramp compensation is performed within a certain time to ensure that the net imbalance power of the zone is continuously transferred to the upper-level plan: in, For partitioning Energy storage Active power command at any given time; The virtual inertia coefficient for energy storage; For partitioning The rate of change of frequency; This refers to the primary frequency regulation ratio for energy storage. For partitioning Frequency deviation; The integral coefficient for energy storage frequency recovery; This represents the maximum available active power capacity for energy storage. This is the sampling step size.

8. A hierarchical distributed control system for grid resilience, employing a hierarchical distributed control method for grid resilience as described in any one of claims 1-7, characterized in that, include: The main site is used to publish scene parameters, global control thresholds, and consistency conditions required for alternating direction multiplier method iteration; Multiple partition edge controllers are configured, with one partition edge controller for each partition. Each partition edge controller has a built-in deterministic real-time kernel, a state estimation module, a local optimization module, and a boundary coordination module. The state estimation module is used to estimate the system state within the partition based on measurement data. The local optimization module is used to solve the local sub-problems corresponding to the partition. The boundary coordination module is used to exchange boundary quantity information with the partition edge controllers of adjacent partitions and perform consistency checks. The measurement unit includes phasor measurement units distributed at various nodes of the power system, used to provide synchronous phasor measurement and time synchronization signals; The execution unit includes an energy storage device, a reactive power compensation device, a controllable load, and a power distribution switchgear, and is used to receive and execute the action commands output by the partition edge controller.

9. A computer-readable storage medium having a program stored thereon, characterized in that, When executed by the processor, the program implements the steps in the hierarchical distributed control method for grid resilience as described in any one of claims 1-7.

10. An electronic device comprising a memory, a processor, and a program stored in the memory and executable on the processor, characterized in that, When the processor executes the program, it implements the steps in the hierarchical distributed control method for grid resilience as described in any one of claims 1-7.

Images