Multi-timescale capacity configuration method for port energy storage systems

Abstract

This invention relates to the field of intelligent port technology and discloses a multi-timescale capacity configuration method for port energy storage systems. The method includes: an offline phase, constructing a safety certificate library containing HEAC segments based on an innovative Hybrid Electromagnetic Arrival Algorithm (HEAC) theory; an online phase, matching early characteristics of synchronous measurement data from grid connection points with the safety certificate library in real time to generate a time-based energy floor curve for real-time control; and finally, applying a Demand Integral Barrier Guarantee (DICB) control law to provide deterministic safety guarantees for supercapacitors, while simultaneously applying the Energy-Deadline Feasibility Release (EDE) criterion for risk-aware energy routing of lithium batteries, and coordinating the regulation of energy storage units based on the energy floor curve to ultimately determine the operational control scheduling table and capacity configuration result list. Thus, this invention achieves safe, economical, and refined configuration and scheduling of port energy storage systems across multiple timescales.

Description

Technical Field

[0001] This invention belongs to the field of port energy storage, and in particular to a method for configuring capacity across multiple time scales in port energy storage systems. Background Technology

[0002] As key hubs for international logistics and energy-intensive areas, the green and electrification transformation of ports has become an irreversible trend. The widespread application of high-power electrical loads such as high-voltage shore power systems, electric container cranes (quay cranes), and electric unmanned trucks in ports has greatly promoted energy conservation and emission reduction.

[0003] However, the integration of these new loads has also brought unprecedented challenges to the port power distribution network. In particular, the millisecond-level inrush currents generated when transformers close, the second-level regenerative energy feedback from quay cranes during operation, and the random start-stop of short-term high-power loads such as ship auxiliary machinery and propellers create complex load characteristics with multiple time scales, high concurrency, and strong impacts. These characteristics easily lead to power quality problems such as voltage dips, flicker, excessive harmonics, and even reverse power flow. In severe cases, this can cause malfunctions in protection devices, threatening the continuity of port operations and the safety and stability of the power grid. Therefore, researching optimized configuration and control methods for energy storage systems suitable for complex port scenarios to mitigate power surges and ensure power quality has significant theoretical and engineering value for supporting the transformation of the port's energy structure and achieving the goal of building a green and smart port.

[0004] Currently, extensive research has been conducted on the capacity configuration and operation control of energy storage systems in power systems. Regarding capacity configuration, existing methods are mostly based on statistical analysis of historical load data from typical days, determining the rated power and capacity of the energy storage system through scenario generation and peak-shaving / valley-filling benefit assessment. For example, opportunity-constrained programming or robust optimization methods are used to optimize the capacity of battery energy storage systems, considering the uncertainties of load and renewable energy forecasts, with the goal of minimizing investment and operating costs. In terms of operation control, research mainly focuses on application scenarios such as energy storage systems participating in grid peak shaving, frequency regulation, and demand-side response. Control strategies are typically based on state feedback, such as setting the charging and discharging power of energy storage according to grid frequency deviation or line power flow, or using model predictive control to continuously optimize the operation strategy over a future period. For hybrid energy storage systems, such as the combination of lithium batteries and supercapacitors, existing energy management strategies typically employ rule-based or filtering methods to decouple high- and low-frequency load power, with the supercapacitor handling the high-frequency components and the lithium battery handling the low-frequency components, leveraging the technological advantages of both.

[0005] However, existing research suffers from profound technical limitations when applied to specific scenarios such as port shore power, where electromagnetic transient shocks are the primary concern. These limitations collectively lead to an imbalance between safety and economy in port energy storage system configuration schemes, making it difficult to achieve truly refined and optimized configurations. Specifically, the core technical problems faced by existing solutions lie in the lack of accurate quantitative assessment methods for the superimposed risks of multi-timescale, asynchronous concurrent electromagnetic transient shocks, and the lack of risk-aware, safety-assured collaborative control mechanisms for hybrid energy storage units. Summary of the Invention

[0006] To address the aforementioned technical problems and to resolve issues such as inaccurate quantification of transient impact superposition risks and inefficient collaborative regulation of hybrid energy storage in existing technologies, this invention provides a multi-timescale capacity configuration method for port energy storage systems.

[0007] The technical solution adopted in this invention is as follows:

[0008] An embodiment of the present invention proposes a multi-timescale capacity configuration method for port energy storage systems, comprising the following steps: establishing an offline safety certificate library, which includes HEAC (Hybrid Electromagnetic Arrival Curve) segments and phase penalty functions generated based on the Hybrid Electromagnetic Arrival Curve theory; matching synchronous measurement data from grid connection points with the safety certificate library online to generate a time-based energy floor curve for real-time control; applying a demand integral barrier-guaranteed control law and an energy-cutoff period feasibility release criterion, and coordinating the regulation of energy storage units based on the time-based energy floor curve to determine the operational control scheduling table and capacity configuration result list.

[0009] The multi-timescale capacity configuration method for port energy storage systems described above in this invention also has the following additional technical features:

[0010] According to one embodiment of the present invention, an offline security certificate library is established, including: performing electromagnetic transient simulations for two events, closing and quay crane feedback, to obtain the power trajectory of the closing segment and the power trajectory of the quay crane feedback segment respectively; extracting the upper envelope of the power trajectory of the closing segment and the power trajectory of the quay crane feedback segment to form the upper bound curve of the closing segment and the upper bound curve of the feedback segment; applying a phase penalty function to define the superposition cost under different phase misalignments, and combining the upper bound curve of the closing segment and the upper bound curve of the feedback segment through a hybrid combination operator to generate a HEAC segment and construct a security certificate library.

[0011] According to one embodiment of the present invention, the phase penalty function is calibrated through the following steps: for a preset set of phase misalignment values, the true mixed arrival upper bound curve of the closing and quay crane feedback events superimposed under the corresponding phase misalignment is obtained; the arithmetic sum of any true mixed arrival upper bound curve and the arrival upper bound curves of the closing segment and the feedback segment under the corresponding phase misalignment is compared to obtain a combined residual sequence; based on the combined residual sequence, the quantile safety margin is used to calculate and fuse measurement and modeling errors to determine the original phase penalty sequence covering the worst case; the original phase penalty sequence is fitted by applying non-negative, symmetric, monotonic, and convex shape constraints to finally constitute the phase penalty function.

[0012] According to one embodiment of the present invention, a demand integral barrier-guaranteed control law is applied to coordinate the regulation of an energy storage unit, including: online matching of grid connection point synchronous measurement data and a safety certificate library, and also used to generate a phase-corrected cumulative demand curve; constructing a demand integral barrier function based on the available electrical energy of the supercapacitor unit, the phase-corrected cumulative demand curve, and a preset upper bound margin of error; solving for a control solution used to maintain the demand integral barrier function in accordance with a preset hold inequality, thereby generating a power reference trajectory for the supercapacitor unit that ensures the energy floor is not crossed.

[0013] According to one embodiment of the present invention, an energy storage unit is coordinated and regulated using an energy-deadline feasibility release criterion, including: online matching of grid connection point synchronous measurement data and a safety certificate library, and further used to discretize the impact event into a set of energy requests with a deadline carrying energy demand and deadline information; for any energy request in the set of energy requests with a deadline, its baseline load ratio is quantified, and a coupling term factor characterizing the risk of the energy request is calculated from three dimensions: slope, energy floor margin, and phase misalignment; based on the combination criterion of the baseline load ratio and the coupling term factor, energy requests that can be safely accepted by the lithium battery energy storage unit are identified, and the remaining energy requests are allocated to the supercapacitor unit to form an energy routing strategy.

[0014] According to an embodiment of the present invention, the process for determining the weight coefficients of the coupling term factors in the combination criterion is as follows: An EDE calibration sample set is constructed, wherein, by simulating and replaying a scenario where any energy request is accepted by a lithium battery energy storage unit, if a hard constraint violation is triggered, the energy request is labeled as unsafe; otherwise, it is labeled as safe. A linear programming problem is established and solved, where the linear programming problem has the hard constraint of blocking all energy requests carrying unsafe labels, and the optimization objective is to minimize the number of rejected energy requests carrying safe labels, thereby obtaining the weight coefficients.

[0015] According to an embodiment of the present invention, the above-mentioned multi-timescale capacity configuration method for port energy storage systems further includes: collecting a set of verification indicators under the execution of the operation period control scheduling table during critical scenario playback or actual operation; comparing the set of verification indicators with the boundary guarantee of the safety certificate library; triggering certificate update when boundary crossing or certificate approximation error exceeds the limit; responding to the certificate update, adjusting the set of parameter uncertainties related to boundary crossing or error exceeding the limit, and re-executing the offline establishment of the safety certificate library to achieve closed-loop tuning of the safety certificate library.

[0016] According to one embodiment of the present invention, the step of online matching of grid-connected point synchronization measurement data and safety certificate library includes: extracting online early feature sequences, including high-frequency voltage distortion, initial magnetization response and closing asymmetry index, from the millisecond-level start segment of grid-connected point synchronization measurement data; performing set consistency matching on HEAC segments in the safety certificate library based on the online early feature sequences to identify the most consistent HEAC segment and estimate the time phase shift; and combining the most consistent HEAC segment and the time phase shift to synthesize a time-type energy floor curve.

[0017] According to one embodiment of the present invention, the combined criterion is specifically manifested as an inequality criterion: the coupling terms of the three dimensions of slope, energy floor margin and phase misalignment are weighted and summed to construct a comprehensive risk coupling term; when the sum of the baseline load ratio and the comprehensive risk coupling term is not greater than a preset release threshold, it is determined that the energy request can be safely accepted by the lithium battery energy storage unit.

[0018] According to an embodiment of the present invention, the steps of determining the operational control scheduling table and the capacity configuration result list include: constructing a joint optimization objective function, which integrates equipment capacity cost, lithium battery aging cost, and voltage and reverse current violation penalties; solving the objective function under the common constraints of the certificate upper bound provided by the safety certificate library, the maintenance inequality of the demand integral barrier guaranteed control law, and the energy-cutoff feasibility release criterion; compiling the optimization result set obtained from the solution into an operational control scheduling table, and deriving a capacity configuration result list containing energy storage unit capacity and operational constraint parameters.

[0019] The beneficial effects of this invention are:

[0020] This invention enables safe, economical, and precise configuration and scheduling of port energy storage systems across multiple time scales. The relevant effects will be described in detail below with reference to specific embodiments. Attached Figure Description

[0021] Figure 1 This is a flowchart of a multi-timescale capacity configuration method for a port energy storage system according to an embodiment of the present invention;

[0022] Figure 2This is a flowchart of offline establishment of a security certificate store according to an embodiment of the present invention;

[0023] Figure 3 This is a flowchart of a calibration phase penalty function according to an embodiment of the present invention;

[0024] Figure 4 This is a flowchart illustrating the application of an energy-deadline feasibility release criterion to collaboratively regulate an energy storage unit according to an embodiment of the present invention.

[0025] Figure 5 This is a flowchart of an energy storage unit that utilizes a demand integral barrier-guaranteed control law for coordinated regulation according to an embodiment of the present invention. Detailed Implementation

[0026] In order to solve the aforementioned problems in the existing technology, the applicant conducted in-depth research and found that:

[0027] Existing methods cannot accurately define the upper bound of the superposition effect of asynchronous concurrent transient events, leading to distorted safety boundaries. Transformer inrush current and quay crane energy feedback are both short-duration electromagnetic transient processes with severe waveform distortion. The peak power and energy resulting from their superposition are closely related to the microsecond to millisecond-level phase misalignments, exhibiting strong nonlinear and non-commutative characteristics. Traditional capacity allocation methods rely on power flow analysis with second- or minute-level resolution, completely failing to capture such transient details. Even using electromagnetic transient simulation to traverse all possible phase misalignments results in exponentially increasing computational costs, rendering it impractical for engineering applications. Simply relying on empirical arithmetic superposition produces overly conservative pseudo-boundaries far exceeding actual physical boundaries, leading to severe over-allocation of energy storage system capacity and wasted investment. This inability to clearly define and explain the risks of transient superposition is the root cause of the dilemma faced by port energy storage planning, which is either unsafe or uneconomical.

[0028] Existing hybrid energy storage control strategies lack the perception and transmission of transient impact risks, failing to achieve true collaborative protection and optimization. Traditional control strategies based on high- and low-frequency decoupling are essentially a form of blind division of labor. They divide power signals solely based on frequency characteristics, without considering the inherent risks—such as the ramp rate of a power request, the overall safety margin of the system at the time of its occurrence, and the nonlinear superposition risks it carries. This one-size-fits-all approach often leads to high-slope, high-risk power surges being incorrectly assigned to lithium batteries, which are highly sensitive to overcurrent and lifespan degradation, accelerating their aging. Simultaneously, it fails to establish a clear and guaranteed energy safety baseline for the supercapacitor's response, making its control objectives vague and unable to provide deterministic safety guarantees. This crude control mechanism renders the collaboration of hybrid energy storage merely a formality, failing to achieve an optimal balance among multiple objectives such as ensuring system safety, protecting equipment lifespan, and reducing operating costs.

[0029] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0030] Example 1: This example describes the overall framework and implementation process of a multi-timescale capacity configuration method for a port energy storage system.

[0031] Specifically, such as Figure 1 As shown, a method for configuring the capacity of a port energy storage system across multiple time scales includes the following steps:

[0032] Step S100: Establish a security certificate library offline. The security certificate library contains HEAC segments and phase penalty functions generated based on the Hybrid Electromagnetic Arrival Calculus (HEMA) theory. Alternatively, establish a security certificate library offline, which contains HEAC segments generated based on the Hybrid Electromagnetic Arrival Calculus (HEMA) theory, and a phase penalty function for combining the HEAC segments.

[0033] Specifically, the safety certificate library is a pre-calculated and stored database. Its core content is not traditional event waveforms or statistical data, but rather an upper bound characterizing the worst-case power demand of a port shore power system under various parameter uncertainties when encountering impact events (such as transformer closing or energy feedback from quay cranes). In this embodiment, this upper bound is generated using an innovative Hybrid Electromagnetic Arrival Algorithm (HEAC) theory. HEAC theory is an improvement and application of arrival curve theory in traditional network calculus, specifically designed to handle the superposition of power impacts in power systems dominated by electromagnetic transient processes, exhibiting strong nonlinearity and phase sensitivity. The establishment of the safety certificate library is completed offline, during the system planning and design phase or the pre-operation preparation phase, utilizing high-performance computing resources for thorough simulation and analysis. This avoids the enormous burden of performing complex electromagnetic transient calculations online in real time, providing a quickly queryable safety guarantee with deterministic boundaries for subsequent real-time control. A HEAC segment is a function segment calculated using this theory, describing the change in worst-case power demand over time after the superposition of one or more types of events. The phase penalty function is a core component of HEAC theory. It is used to quantify the additional cost of the superposition effect caused by different impact events occurring at different times (i.e., phase misalignment). Its specific composition will be detailed in subsequent embodiments.

[0034] Step S200: Online matching of grid-connected point synchronous measurement data with the safety certificate library to generate a time-based energy floor curve for real-time control. Alternatively, online extraction of early feature sequences from the grid-connected point synchronous measurement data, followed by matching and phase alignment using the early feature sequences in the safety certificate library to generate a time-based energy floor curve and a phase-corrected cumulative demand curve.

[0035] In this embodiment, when the port shore power system is actually running, it collects synchronous measurement data such as voltage and current at the grid connection point in real time through a high-precision, high-sampling-rate synchronous phasor measurement unit (PMU) or similar device. The online matching step involves performing a rapid pattern recognition and matching between the early features of these millisecond-level real-time measurement data and HEAC segments in an offline security certificate library. It is understood that different initial conditions (such as transformer residual flux, upstream grid short-circuit capacity, etc.) will lead to differences in the waveform characteristics of the impact event. By matching the early features of the online data, the HEAC segment most consistent with the current scenario can be quickly identified from the library, and the precise time phase of the event can be estimated. Based on this, a time-based energy floor curve is generated. This curve is not a fixed energy threshold, but a curve that dynamically changes over time. It defines the minimum available energy level that the energy storage system must maintain to ensure system safety and stability and prevent voltage overruns or reverse power flow problems within a future period (e.g., from milliseconds to seconds) starting from the current moment. This curve is the core basis for subsequent coordinated control of the energy storage unit.

[0036] Step S300: Applying the demand integral barrier guarantee control law and the energy-deadline feasibility release criterion, and coordinating the regulation of energy storage units based on the time-based energy floor curve, thereby determining the operational control scheduling table and capacity configuration result list. Alternatively, based on the time-based energy floor curve and the phase-corrected cumulative demand curve, and in conjunction with the safety certificate library, applying the demand integral barrier guarantee control law and the energy-deadline feasibility release criterion, coordinating the power allocation between the supercapacitor unit and the lithium battery energy storage unit, and generating the operational control scheduling table.

[0037] After obtaining the time-based energy floor curve, this step aims to differentiate and coordinate the regulation of different energy storage units (exemplarily including supercapacitor energy storage units and lithium battery energy storage units) in the port energy storage system. It is understood that supercapacitors have high power density and fast response, making them suitable for handling millisecond-level severe power surges, but their energy capacity is limited and their cost is high; while lithium batteries have high energy density and relatively low cost, but they are more sensitive to high-rate, high-frequency charging and discharging, which accelerates their aging. Therefore, this embodiment employs two dedicated control strategies: First, a Demand Integral Control Barrier (DICB) guaranteed control law applied to the supercapacitor energy storage unit. This control law treats the time-based energy floor curve as an insurmountable barrier, actively adjusting the power output of the supercapacitor through a guaranteed control method to ensure that its available energy remains above the floor curve at any given time, thereby providing deterministic safety assurance for the system. Second, an Energy Deadline Feasibility Emission (EDE) criterion applied to the lithium battery energy storage unit. This criterion decomposes complex impact power into a series of micro-tasks with energy demands and deadlines, and comprehensively assesses the risk of lithium batteries undertaking these tasks from multiple dimensions, such as slope, energy floor margin, and phase misalignment risk. Only energy requests deemed safe are allowed to be processed by lithium batteries; other requests with higher risks or faster response requirements are forced to be handled by supercapacitors. Through the synergistic effect of the DICB control law and the EDE criterion, the system can maximize the use of lithium batteries to handle low-risk, slow-dynamic energy demands while ensuring safety, reducing the configuration capacity and number of supercapacitor operations, and protecting lithium batteries from high-frequency impacts, thus reducing their aging costs. Finally, through optimization within an operating cycle (e.g., a day or a ship berthing cycle), considering factors such as equipment capacity cost, aging cost, and power quality violation penalties, the optimal operating period control schedule (i.e., power commands and parameter settings for each energy storage unit) is determined, and a capacity configuration result list containing the final capacity selection and operating constraint parameters of the energy storage units is derived.

[0038] According to one aspect of this application, it further includes: replaying and verifying the execution effect of the runtime control scheduling table, and updating the security certificate database when there is a deviation between the execution indicator and the upper bound of the security certificate database.

[0039] To enable the above method, one specific implementation of this embodiment also includes a process of aggregating and modeling basic data. For example, by aggregating and aligning grid-connected point synchronous measurement data, equipment and parameter ledgers, access and switching event records, operation and maintenance and protection records, technical limits and compliance constraints, etc., the timestamps and phases are unified to form a synchronization data packet. Based on this, based on the synchronization data packet, time windows are segmented according to scenarios, such as millisecond-level time windows (covering the instant of point wave position closing), second-level time windows (covering excitation inrush current attenuation), and ten- to thirty-second time windows (covering short-term power fluctuations of the thruster and quay crane), to obtain time window indices. Further, within each window, closing impact power trajectories, excitation inrush envelope characteristics, quay crane feedback power curves, etc., are generated from the grid-connected point synchronous measurement data, and aggregated into a windowed dataset, providing high-quality data input for subsequent offline modeling and online matching.

[0040] Example 2: This example elaborates on the offline method for establishing the S100HEAC security certificate library in Example 1.

[0041] In this embodiment, as Figure 2 As shown, the steps for establishing a security certificate store offline include:

[0042] Step S110: Perform electromagnetic transient simulation for the two events of closing and quay crane feedback, and obtain the power trajectory of the closing segment and the power trajectory of the quay crane feedback segment respectively.

[0043] Specifically, this step involves reproducing high-frequency impact events in the port shore power system by establishing a high-fidelity electromagnetic transient simulation model. The closing event specifically refers to the closing process of the shore power system transformer supplying power to berthing vessels, a process that generates intense inrush current. The quay crane feedback event specifically refers to the process where the motor of the quay crane is in generator mode when lowering heavy objects such as containers, feeding regenerative energy back to the grid. These two events are the most significant sources of impact load in the port shore power system. The simulation model will cover the nonlinear magnetization and hysteresis characteristics of the transformer, the equivalent impedance of the upstream grid, the distributed parameters of the connecting cables, and key equipment such as the point-wave closing device and pre-insertion resistors. Through extensive, batch-processed sweep simulations within a hyperparameter space containing various uncertain parameters (such as transformer closing angle, residual magnetization level, upstream grid short-circuit capacity, and the proportion of quay crane feedback power), a series of original power time series representing different operating conditions can be obtained, namely, the power trajectory of the closing segment and the power trajectory of the quay crane feedback segment.

[0044] Step S120: Extract the upper envelope of the power trajectory of the closing segment and the power trajectory of the quay crane feedback segment to form the upper boundary curve of the closing segment and the upper boundary curve of the feedback segment.

[0045] After obtaining the original power trajectory, an upper envelope extraction is required to form an arrival curve with a deterministic upper bound. In this embodiment, a non-causal upper envelope extraction algorithm combining sliding maximum and convex hull is preferably used. Sliding maximum refers to taking the maximum power within a sliding time window to form a preliminary envelope; convex hull combination corrects the preliminary envelope to ensure that the final generated arrival upper bound curve is convex upwards. This property is crucial for the convergence and simplification of subsequent calculations. Since this step is performed offline, non-causal processing can be used. That is, when extracting the envelope value at time t, data after time t can be used to obtain a more compact and accurate upper bound than online causal filtering. Through this step, a series of discrete, fluctuating power trajectories are transformed into a set of function curves with a clear mathematical form (e.g., piecewise linear or piecewise affine) that guarantee that the actual power will not exceed the limit under any circumstances, namely, the arrival upper bound curve of the closing segment and the arrival upper bound curve of the feedback segment.

[0046] Step S130: Apply the phase penalty function to define the superposition cost under different phase misalignments, and combine the upper bound curve of the closing segment and the upper bound curve of the feedback segment by the hybrid combination operator to generate HEAC segment and build a security certificate library.

[0047] It is understandable that when a closing event and a quay crane feedback event occur simultaneously or successively, the total power impact is not a simple arithmetic superposition, but is closely related to the relative difference in the timing of the two events (i.e., phase misalignment). To accurately characterize this complex superposition effect, this embodiment introduces a hybrid combination operator and a phase penalty function.

[0048] The specific application of the hybrid combination operator is defined as follows: The infimum is calculated over the range of all possible phase misalignments to determine the HEAC segment; the value corresponding to the infimum is equal to the minimum sum of three terms: the closing segment reaching the upper bound curve, the feedback segment reaching the upper bound curve after time shifting with the corresponding phase misalignment, and the phase penalty function under that phase misalignment. Its mathematical expression can be written as: A mix (t)=inf Δ∈R (A close (t)+A regen (t-Δ)+ (Δ));

[0049] Among them: A mix (t) represents the final generated HEAC segment, indicating the worst-case power upper bound at time t after the superposition of the two events, in watts (W); t is the time variable, in seconds (s); A close The curve for the closing segment reaching the upper bound is a function of time; A regen(t-Δ) represents the curve of the feedback segment reaching the upper bound after time shift; Δ is the phase misalignment variable, representing the time delay of the feedback event relative to the closing event, in seconds (s). (Δ) is the phase penalty function, which is a function of phase misalignment Δ. It quantifies the additional superposition power caused by the phase misalignment, and its unit is watts (W). Its specific calibration method is detailed in another embodiment; inf Δ∈R This represents the infimum of the expression over all possible real phase misalignments Δ. Through this calculation, a unified upper bound for superposition events, i.e., the HEAC segment, can be obtained that no longer depends on the specific phase misalignment Δ.

[0050] Furthermore, the hybrid combination operator also satisfies the subadditivity property: for an event sequence consisting of two consecutive events, the overall mixture reaches an upper bound, which is no greater than the sum of the individual mixtures of each event in the sequence reaching their upper bounds, minus a phase alignment gain resulting from sequence concatenation. This property A mix 1⊕2 (t)≤A mix 1 (t)+ A mix 1 (t)(t)-ψ(Δ12), where ψ(Δ12) is the phase alignment gain, makes it easy to extend this theory to the analysis of multiple event sequences.

[0051] The HEAC segments generated for different initial operating conditions and their corresponding metadata (such as applicable parameter ranges, error upper bounds, etc.) are stored together in the database to form a safety certificate library.

[0052] Furthermore, this embodiment also includes: performing time integration on the HEAC segment in the security certificate repository to obtain a cumulative demand curve characterizing the cumulative energy demand under the worst-case scenario. That is, E cum (t)= 0 t mix (τ)dτ;where E cum (t) is the cumulative demand curve. mix (τ) represents the HEAC segment.

[0053] Optionally, based on this, the cumulative demand curve can be sampled in a scenario by combining the set of parameter uncertainties, and quantiles can be taken according to a preset confidence level (e.g., 95% or 99%) to generate a quantile energy floor baseline curve, providing a benchmark for generating a time-type energy floor in the online phase.

[0054] Example 3: This example is a further refinement of Example 2, elaborating on the phase penalty function. The specific calibration and construction process of (Δ). It is understandable that the accuracy of this function is directly related to the compactness and security of HEAC fragments.

[0055] Specifically, such as Figure 3 As shown, the phase penalty function is calibrated through the following steps:

[0056] Step S210: For the preset set of phase misalignment values, obtain the true mixed curve of the closing and quay crane feedback events superimposed under the corresponding phase misalignment to reach the upper limit.

[0057] Specifically, this step aims to construct a high-quality ground truth dataset for function fitting. First, a discrete set of phase misalignment values ​​covering the typical working range needs to be predefined; exemplarily, this set can be defined as Δ. grid ={0, ±1ms, ±2ms, ..., ±60ms}, where the maximum phase misalignment can be selected between 40ms and 100ms according to the actual operating cycle of the port, and the step size can be adjusted between 0.5ms and 2ms. Subsequently, for each phase misalignment value Δ_k in this set, and for each hyperparameter grid point in Example 2 (i.e., different initial operating condition combinations), a complete superimposed event simulation is performed using a high-fidelity electromagnetic transient model. In this simulation, the closing event occurs at t=0, while the quay crane feedback event occurs precisely at t=Δk. By extracting the upper envelope of the total grid connection point power trajectory obtained from the simulation (for example, using the same sliding maximum and convex hull combination algorithm as in Example 2), the true mixed arrival upper bound curve Atrue under this specific operating condition and specific phase misalignment can be obtained. _mix (t;Δk). Simultaneously, the upper bound curve A for the case with only a closing event was also independently recorded in this simulation. close (t) and the upper bound curve A when there are only feedback events. regen (t). These three curves, along with their corresponding phase misalignment values ​​Δk and operating condition metadata, are used as a sample. This process is repeated until all preset operating conditions and phase misalignment values ​​are traversed, ultimately forming a phase misalignment sample library.

[0058] Step S220: Compare any real mixed arrival upper bound curve with the arithmetic sum of the arrival upper bound curves of the closing segment and the feedback segment under the corresponding phase misalignment, and quantize to obtain the combined residual sequence.

[0059] The purpose of this step is to quantify the difference between the true superposition effect and a simple arithmetic sum; this difference is the nonlinear superposition cost caused by phase misalignment. For each sample in the sample library, at each discrete time point t, its combined residual is calculated using the formula: r(t;Δ)=Atrue mix(t;Δ)-(A close (t)+A regen (t-Δ)); where: r(t; Δ) is the combined residual at time t with a phase misalignment of Δ, in watts (W); A true_mix (t; Δ) represents the upper bound curve of the true mixture for this sample; A close (t) represents the curve of the corresponding closing segment reaching the upper bound; A regen (t-Δ) represents the curve after shifting the corresponding feedback segment to the upper bound curve along the time axis by Δ. It can be seen that the physical meaning of r(t; Δ) is the minimum amount of power that needs to be compensated at time t in order to encompass the true superposition power using simple arithmetic. Calculating this for each sample over the entire time domain yields a time-varying sequence of combined residuals.

[0060] Step S230: Based on the combined residual sequence, the phase penalty original sequence covering the worst case is determined by calculating and fusing measurement and modeling errors using quantile safety margin.

[0061] After obtaining the residual sequence for each sample, a conservative upper bound that can cover all operating conditions and all time points needs to be extracted. First, for a given phase misalignment value Δk, all corresponding samples are found from the sample library. For each sample's residual sequence r(t; Δk), the supremum (i.e., peak value) is taken at time t, resulting in a set of residual peak values. Then, a quantile safety margin is calculated on this set of peak values, for example, taking its 99th quantile, thus obtaining a residual upper bound r that statistically covers the vast majority of worst-case scenarios. (Δk)=Q 0.99 max t r(t; Δk)); where Q 0.99 99 raw (Δk)=r (Δk) + δmeas + δmodel; where: raw (Δk) represents the original phase penalty value corresponding to the phase misalignment Δk, which is a point in the discrete sequence; δmeas is the equivalent power compensation margin calculated from sensor measurement delay, data synchronization error, etc.; δmodel is the power compensation calculated from the upper bound of the modeling error caused by the difference between the electromagnetic transient model and the actual physical system. Repeating this calculation for each phase misalignment value in Δgrid yields a discrete original phase penalty sequence covering the worst-case scenario. raw (Δ).

[0062] Step S240: Apply non-negative, symmetric, monotonic, and convex shape constraints to the original phase penalty sequence for fitting processing, and finally construct the phase penalty function.

[0063] The original phase-penalized sequence is discrete and may contain noise or fluctuations that do not satisfy certain physical properties. To obtain a continuous, smooth, and well-behaved function for subsequent calculations, it needs to be fitted. In this embodiment, the fitting process is constructed as an optimization problem with a set of shape constraints. The shape constraints include:

[0064] Nonnegativity (Δ)≥0, and (0) = 0, because there should be no penalty when there is no phase misalignment.

[0065] symmetry, (Δ)= (-Δ), because the impact of a feedback event being earlier or later by the same amount of time should be symmetrical.

[0066] Monotonicity, when |Δi|<|Δj|, should have (Δi)≤ (Δj), that is, the greater the phase misalignment, the greater the uncertainty and superposition risk, and the penalty function should be monotonically undecreasing.

[0067] Convexity ensures that the second derivative of the function is non-negative. This constraint guarantees that the function has only one global minimum and makes the optimization problem of finding the infimum in step S130 of Example 2 a convex problem, guaranteeing the uniqueness and efficiency of the solution. Under these constraints, the goal is to minimize the fitting error (e.g., ∑ k ( (Δk)- raw (Δk)) 2 With the objective function , an optimal piecewise convex function can be found. (Δ) (e.g., a piecewise linear or piecewise quadratic function). This function is the final calibrated phase penalty function that can be written into the security certificate store. For example, a calibrated function might, after normalization, have a value at the maximum phase misalignment Δmax. norm (Δ max It falls within the range of 0.05 to 0.30.

[0068] Example 4: This example elaborates on step S200 of Example 1, the online matching and time-type energy floor generation method, and describes how to apply the offline-constructed theoretical model to online real-time control.

[0069] In this embodiment, the step of online matching of grid connection point synchronization measurement data and security certificate library includes:

[0070] Step S310: Extract the online early feature sequence, including high-frequency voltage distortion, initial magnetization response and closing asymmetry index, from the millisecond-level starting segment of the synchronous measurement data at the grid connection point.

[0071] When the system detects a suspected impact event (e.g., a circuit breaker status change signal or a sudden current surge), it immediately initiates a millisecond-level time window for high-sampling-rate data acquisition. The purpose of this step is to quickly diagnose the inherent properties of the event within a very short time after its occurrence (e.g., the first 1-3 power frequency cycles, approximately 20-60ms). It is understood that key parameters such as residual magnetic flux inside the transformer and the instantaneous short-circuit capacity of the upstream power grid cannot be directly measured in real time, but they are reflected like fingerprints in the details of the voltage and current waveforms in the early stages of the event. The online early feature sequence extracted in this embodiment is a quantitative representation of these fingerprints. Specifically, the feature sequence may include:

[0072] High-frequency voltage distortion index: This index is calculated by performing a short-time Fourier transform (STFT) or wavelet transform on the grid connection point voltage signal, determining the energy percentage or peak amplitude within the high-frequency band from 1kHz to 2kHz. This index is highly sensitive to the saturation level of the transformer core.

[0073] Initial magnetization response index: This index calculates the DC component offset of the excitation current and the asymmetry of the three-phase current peak values ​​during the first few cycles after closing the circuit. This index directly reflects the residual magnetic flux level and phase angle at the instant of closing the circuit.

[0074] Closing asymmetry index: This index calculates the deviation of the phase difference between the three-phase voltage or current and the ideal 120 degrees, as well as the amplitude of the instantaneous zero-sequence component. This index reflects the synchronicity of three-phase closing and the background imbalance of the system. Combining these calculated indices into a feature vector constitutes the online early characteristic sequence.

[0075] Step S320: Based on the online early feature sequence, perform set consistency matching on HEAC fragments in the security certificate library to identify the most consistent HEAC fragment and estimate the time phase shift.

[0076] This step serves as a bridge between theory and reality. First, the online early feature sequence vector extracted in step S310 is compared with the metadata (i.e., the theoretical feature vector corresponding to the operating parameters used to generate the segment) attached to each HEAC segment in the security certificate repository. This comparison can be a set consistency matching, for example, calculating the Euclidean or Mahalanobis distance between the online feature vector and each theoretical feature vector in the repository. The HEAC segment corresponding to the smallest distance is identified as the HEAC segment most consistent with the current actual scenario. While identifying the most consistent segment, it is also necessary to accurately estimate the start time of the event. This can be achieved by performing cross-correlation or least-squares fitting between the actually measured current waveform (e.g., the first 100ms) and the theoretical typical waveform corresponding to the most consistent HEAC segment. The optimal time offset estimate Δτ can be obtained by solving Δτ = arg min. δ ||Φmeas(t) - ΦHEAC(t-δ)||2 is obtained; where: Δτ is the estimated time phase offset; Δ is the candidate time offset used for searching; Φmeas(t) is the actual measured waveform sequence; ΦHEAC(t-δ) is the theoretical waveform sequence corresponding to the most consistent HEAC segment;|| ||2 represents the L2 norm (i.e., Euclidean distance).

[0077] In some alternative implementations, this matching process can also be implemented using a pre-trained machine learning model (such as a K-nearest neighbor classifier, support vector machine, or small neural network) that takes the early feature sequence as input and directly outputs the index number and temporal phase shift of the most consistent HE-AC segment to achieve a faster matching speed.

[0078] Step S330: Combine the most consistent HEAC segment with the time phase shift to synthesize a time-type energy floor curve.

[0079] After determining the most consistent HEAC segment and the precise time phase offset Δτ, a forward-looking safety boundary can be tailored for the current event. The synthesis process specifically involves: using the power upper bound function A represented by the selected most consistent HEAC segment... mix (t), based on the estimated phase shift Δτ, a power upper bound prediction curve A aligned with the current event is obtained by time shifting. predicted (t)= A mix (t-Δτ). Then, for this predicted curve from the current time t... now By integrating forward over time, we can obtain the cumulative energy demand forecast curve E. demand (t)= tnow t A predicted(τ)dτ. This E demand The (t) curve is the final time-based energy floor curve used for real-time control. It dynamically depicts the minimum energy value that the energy storage system needs to reserve at each future moment in order to fully absorb the upcoming power surge, providing a precise and reliable target trajectory for subsequent guaranteed control.

[0080] Example 5: This example elaborates on one of the collaborative control strategies in step S300 of Example 1. The Energy-Deadline Feasibility Release Criterion (EDE) is explained in detail. A refined decision-making mechanism is established to identify energy requests with different risk characteristics, thereby achieving effective protection and efficient utilization of the lithium battery energy storage unit.

[0081] In this embodiment, as Figure 4 As shown, the steps for coordinating the regulation of energy storage units using the energy-deadline feasibility release criterion include:

[0082] Step S410 involves online matching of grid connection point synchronous measurement data and safety certificate library, and also discretizing the impact event into a set of deadline energy requests carrying energy demand and deadline information.

[0083] Specifically, in Example 4, a power upper bound prediction curve A aligned with the current event was generated. predicted After (t), to facilitate discretized decision-making and scheduling, this continuous power curve needs to be transformed into a set of independent, manageable task units. This discretization process can be performed by adjusting A along the time axis. predicted (t) This is achieved by micro-window slicing; for example, the time window width of the slice can be set to 1 millisecond to 5 milliseconds. For the i-th time window [t_i, t_{i+1}], its corresponding energy demand is calculated as the time integral of the power curve within that window, i.e., energy. i = [t i ,t i+1 A predicted (τ)dτ. The end time t of this window. i+1 This is defined as the deadline_i for the energy demand. By slicing the power prediction curve of the entire event, a continuous impact event is transformed into a set of multiple {energy_i, deadline_i} tuples with a time-sequential relationship, namely the deadline energy request set. This transformation creatively maps a transient power absorption problem in a power system to a task scheduling problem in the field of real-time computing, laying the foundation for the introduction of more advanced scheduling and decision-making algorithms.

[0084] Step S420: For any energy request in the energy request set at the deadline, quantify its baseline load ratio, and calculate the coupling term factor characterizing the risk of the energy request from three dimensions: slope, energy floor margin, and phase misalignment.

[0085] For any request i in the energy request set with a deadline, to comprehensively assess its suitability for handling by lithium batteries, it needs to be quantitatively characterized from multiple dimensions. First, its baseline load ratio is quantified. This metric primarily assesses the ratio of the requested energy to the battery's processing capacity. It is calculated as: baseline_load_ratio_i = energy_i / available_release_i; where: baseline_load_ratio_i is the baseline load ratio of the i-th request, a dimensionless value; energy_i is the energy demand of the request, in joules (J); available_release_i is the maximum energy that the lithium battery cell can safely release based on its current state of health (SOH), state of charge (SOC), temperature, and other constraints up to its deadline_i. Second, a coupling term factor characterizing the risk of this energy request is calculated. This factor is not a single indicator but a comprehensive assessment term integrating information from three key risk dimensions:

[0086] Slope: Quantifies the intensity of the power surge request. Specifically, it can be calculated as the rate of change or ratio of the average power (energy_i / (t_{i+1}-t_i)) within the time window containing request i relative to the maximum safe charge / discharge power allowed by the lithium battery. A steep power rise (high slope) is extremely detrimental to the electrochemical system and lifespan of a lithium battery.

[0087] Energy floor margin: Quantifies the safety redundancy of the entire energy storage system when a request occurs. Specifically, it can be calculated as the minimum difference between the actual available energy held by the supercapacitor unit within the time window when request i occurs and the time-based energy floor curve generated in Example 4. If this margin is small, it indicates that the system is near the safety boundary, and any additional disturbance (such as a lithium battery response that is not as expected) may cause system instability, increasing the risk of having the lithium battery take over the task.

[0088] Phase misalignment: An inherent nonlinear risk inherent in the quantization request itself. This risk is inherited from the HEAC segment that generated the request, and its value can be directly taken from the phase penalty function calibrated in Example 3. \phi(\Delta) The function value under phase misalignment of the current event. A large phase penalty value means that the superposition effect of the current event is very complex and severe.

[0089] Step S430: Based on the combined criteria of baseline load ratio and coupling term factor, identify energy requests that can be safely handled by the lithium battery energy storage unit, and allocate the remaining energy requests to the supercapacitor unit to form an energy routing strategy.

[0090] Specifically, the combined criterion is manifested as an inequality criterion. First, the risk factors calculated in step S420—slope, energy floor margin, and phase misalignment—are weighted and summed to construct a comprehensive risk coupling term Γ. Then, the baseline load ratio is summed with this comprehensive risk coupling term and compared with a preset release threshold (e.g., threshold 1). When the sum of the baseline load ratio and the comprehensive risk coupling term is not greater than the preset release threshold, it is determined that the energy request can be safely accepted by the lithium battery energy storage unit. The complete form of this inequality criterion can be written as:

[0091] baseline_load_ratio_i + Γ(slope_ratio_i, floor_deficit_i, phase_penalty_i) ≤ 1; where: Γ(·) is the comprehensive risk coupling term function, which in this embodiment is preferably in linear weighted form, i.e.: Γ = a1 slope_ratio_i+a2 floor_deficit_i+a3 `phase_penalty_i`; where: a1, a2, a3 are three non-negative weight coefficients, corresponding to the importance of the three risk dimensions of slope, floor margin, and phase misalignment, respectively, and their specific calibration methods are detailed in the next embodiment; `slope_ratio_i`, `floor_deficit_i`, and `phase_penalty_i` are the quantified values ​​of the three risk dimensions after normalization. For each request in the request set, this criterion is calculated once. If the inequality is true, a scheduling instruction is generated to route the energy request `energy_i` to the controller of the lithium battery energy storage unit; if the inequality is false, the request is determined to be a high-risk request and is forcibly routed to the supercapacitor energy storage unit for absorption. By judging the order of all requests, a complete and differentiated energy routing strategy is finally formed.

[0092] In some alternative implementations, the comprehensive risk coupling term function Γ(·) can also take a nonlinear form, for example, through a small artificial neural network or a multidimensional lookup table generated from a large amount of offline simulation data, to capture more complex nonlinear coupling relationships between risk dimensions. Furthermore, the weighting coefficients a1, a2, a3 can also be dynamic coefficients that are adaptively adjusted according to the real-time updated state of health (SOH) of the lithium battery. For example, when the battery is severely aged, the weight a1 of the slope risk is increased accordingly.

[0093] Example 6: This example, as a further refinement and support of Example 5, describes the calibration method for the weights of the coupling terms in the EDE criterion, especially how the core parameters in the combined criterion, the weight coefficients a1, a2, and a3, are calibrated offline using a systematic method that prioritizes safety.

[0094] The weighting coefficients of the coupling term factors in the combined criterion are determined through the following steps:

[0095] Step S510: Construct an EDE calibration sample set. In this set, a scenario in which any energy request is handled by a lithium battery energy storage unit is replayed through simulation. If a hard constraint violation is triggered, the energy request is labeled as unsafe; otherwise, it is labeled as safe.

[0096] This step aims to generate truth-labeled training data for subsequent supervised learning through simulation experiments. First, a large set of deadline energy requests is collected or generated from historical data or large-scale simulations. For each energy request *i*, a scenario replay simulation assuming the lithium battery will handle the request is performed. In this simulation, the states of other parts of the system (such as the power grid, supercapacitors, etc.) remain consistent with when the request occurs, but the requested energy *energy_i* is forcibly injected as an instruction into the simulation model of the lithium battery energy storage unit. Throughout the simulation, several key system-level hard constraints are closely monitored, such as: whether the voltage deviation at the grid connection point exceeds the compliance range (e.g., ±5%), whether power reverse current occurs, and whether the internal state of the lithium battery itself triggers protection thresholds (e.g., cell current rate, temperature, terminal voltage, etc.). If any hard constraint is triggered or violated during the simulation, it is considered unsafe for the lithium battery to handle the request, and an unsafe label is assigned to request *i* (e.g., label_i=1). Conversely, if the simulation completes successfully and all constraints remain within safe limits, a safety label (label_i=0) is assigned to the request. Simultaneously, for each labeled request, its corresponding feature vector {baseline_load_ratio_i, slope_ratio_i, floor_deficit_i, phase_penalty_i} is calculated and stored. This process is repeated to ultimately form a richly illustrated EDE calibration sample set containing a large number of samples.

[0097] Step S520: Establish and solve a linear programming problem. The linear programming problem has the hard constraint that all energy requests carrying unsafe tags must be blocked, and the optimization objective is to minimize the number of rejected energy requests carrying safe tags, thereby obtaining the weight coefficients.

[0098] After obtaining the calibration sample set, this step aims to solve for an optimal set of weight coefficients {a1, a2, a3}. Here, "optimal" does not refer to the highest classification accuracy in the traditional sense, but rather prioritizes safety as the absolute highest priority. Therefore, the problem is constructed as a special linear programming (LP) problem. The decision variables for this linear programming problem are the weight coefficients a1, a2, a3. Its constraints include:

[0099] Hard constraint on unsafe samples (zero miss constraint): For all samples j with a label of 1 (unsafe), the result of their criterion calculation must be forced to be greater than the allowance threshold of 1, i.e.: L_j + a1x1_j + a2x2_j + a3 x3_j≥1+ε. Where: L_j, x1_j, x2_j, x3_j are the baseline load ratio of sample j and three normalized risk factors, respectively; ε is a small positive safety margin (e.g., 0.05) used to increase the robustness of the criterion. This constraint ensures that all known dangerous operations will be prevented by the criterion.

[0100] Soft constraint for safe samples: For all samples k with a label of 0 (safe), the calculated criterion result should be as small as possible less than the release threshold of 1, i.e.: L_k + a1x1_k + a2x2_k + a3 x3_k≤1-ε+ξ_k. Where: ξ_k(xi_k) is a non-negative slack variable.

[0101] Non-negativity constraints for coefficients: a1≥0, a2≥0, a3≥0. This aligns with physical intuition, meaning that increasing risk factors only increases total risk, not decreases it.

[0102] The optimization objective of this linear programming problem is to minimize the sum of slack variables for all safe samples, i.e., minΣ_ k ξ_ k Minimizing the sum of slack variables is equivalent to ensuring that the criterion result for as many safe samples as possible is less than the threshold 1, while satisfying all the hard constraints mentioned above. In other words, it minimizes the number of operations that could have been safely executed but were incorrectly rejected. By solving this linear programming problem, we can obtain a set of optimal, essentially safety-prioritized, weight coefficients a1, a2, a3.

[0103] In some preferred embodiments, to prevent the criterion from fluctuating near the threshold, a hysteresis band can be introduced into the criterion. For example, the state is switched only when the previous decision passed and the current value is <1-width_hys, or the previous decision failed and the current value is >1+width_hys, where width_hys is the hysteresis width (e.g., 0.02). Furthermore, the obtained weight coefficients a1, a2, a3, and related thresholds, normalization parameters, etc., will be written as part of the metadata into the security certificate store for use by the online system.

[0104] Example 7: This example elaborates on one of the collaborative control strategies in step S300 of Example 1, the Demand Integral Barrier Guarantee Control Law (DICB). A control strategy is designed for the supercapacitor unit with safety as the primary objective, ensuring that the available energy of the energy storage system will not fall below the dynamically changing safety baseline under any circumstances, thus providing a solid guarantee for the stable operation of the entire port shore power system.

[0105] In this embodiment, as Figure 5 As shown, the steps for coordinating the control of an energy storage unit using a demand integral barrier-guaranteed control law include:

[0106] Step S610 involves online matching of grid connection point synchronization measurement data with the security certificate library, and also generates a phase-corrected cumulative demand curve.

[0107] Specifically, in Example 4, after identifying the most consistent HEAC segment through online matching and estimating the time phase shift, a power upper bound prediction curve A that is precisely aligned with the current impact event can be obtained. predicted (t). By examining this curve from the current time t... now By performing time integration, we obtain the phase-corrected cumulative demand curve, denoted as E. cum (t). The physical meaning of this curve is: in order to completely absorb and suppress the impact event currently occurring, from the current moment to any future moment t, the cumulative energy demand that the system is expected to face under the worst-case scenario. This curve constitutes a dynamic and forward-looking energy impact trajectory, which is the defense target that the DICB control law needs to address.

[0108] Step S620: Based on the available electrical energy of the supercapacitor unit, the cumulative demand curve after phase correction, and the preset upper bound of error, construct the demand integral barrier function.

[0109] This step is the core idea of ​​the DICB method. Traditional control methods typically track a given power or voltage reference value, while the DICB control law defines a safe operating domain for the system by constructing a barrier function, ensuring that the system's state trajectory never crosses the boundary of this safe domain. In this embodiment, the demand integral barrier function is defined as: B(t) = E_supercap(t) - Ecum(t) - M; where: B(t) is the barrier function value at time t, its physical meaning is the safety margin of the current supercapacitor's energy reserve relative to the worst-case energy demand, in joules (J), and the system is in a safe state when B(t) > 0; E_supercap(t) is the actual available electrical energy of the supercapacitor unit at time t, which can be determined by real-time measurement of its terminal voltage V_sc(t) and the formula E_supercap(t) = 0.5. C_sc The value is calculated as (V_sc(t)^2-V_min^2), where C_sc is the total capacitance of the supercapacitor and V_min is the minimum allowable operating voltage of the supercapacitor; E cum(t) is the cumulative demand curve after phase correction generated in step S610; M is a preset, always positive upper bound error margin in joules (J). It is set to cope with various uncertainties, including the approximation error of the HEAC model itself, the delay of sensor measurement and data transmission, and the lag in the execution of control commands, so as to make the safety boundary more conservative and robust.

[0110] Step S630: Solve for the control solution used to maintain the demand integral barrier function in accordance with the preset hold inequality, thereby generating a power reference trajectory for the supercapacitor unit that ensures the energy floor is not crossed.

[0111] After constructing the barrier function B(t), the control objective becomes: to maintain B(t) ≥ 0 throughout the entire event duration. To achieve this objective, the DICB control law imposes a pre-defined hold inequality. In the continuous-time domain, this inequality can be expressed as: (t)≥-k B(t); where: (t) represents the first derivative of the barrier function B(t) with respect to time, i.e., the rate of change of the safety margin; k is a positive constant, called the feedback coefficient or convergence rate. This inequality means that the safety margin B(t) is allowed to decrease, but its rate of decrease cannot exceed k times its own size. In other words, when the system state is far from the safety boundary (B(t) is large), the margin is allowed to decrease at a relatively fast rate; while when the system state approaches the safety boundary (B(t) tends to 0), the rate of decrease of the margin is also forced to tend to 0, thus forming a strong repulsive force at the boundary, preventing the state trajectory from crossing the boundary. In a digital control system discretized with a control period T_s, the above inequality is equivalent to: B[k+1]-(1-k T_s) B[k]≥0; where B[k] and B[k+1] are the barrier function values ​​for the current control cycle and the next control cycle, respectively. Since the value of B[k+1] is directly related to the supercapacitor power P_supercap[k] in the next cycle, this inequality constitutes a core constraint on P_supercap[k]. In each control cycle k, the controller solves a small-scale, real-time optimization problem to determine the optimal power reference value P. _supercap [k]. This optimization problem can be defined as minimizing the change in power command (i.e., min(P)) while satisfying the discrete-preservation inequality and the physical constraints such as the power ramp-up rate and voltage operating range of the supercapacitor itself. _supercap [k]-P _supercap [k-1]) 2This ensures the smoothness of control. By continuously solving this optimization problem, a series of power reference values ​​are generated, which are the power reference trajectories that ensure the energy floor is not crossed.

[0112] In some alternative implementations, the feedback coefficient k can be an adaptive parameter. For example, the value of k can be dynamically increased when a large system disturbance or an increase in model prediction error is detected, to provide stronger safety assurance. Furthermore, the barrier function can also take other forms, such as a logarithmic barrier function or an exponential barrier function, to obtain different control characteristics when approaching the boundary.

[0113] Example 8: This example further supplements and optimizes the dual-slope release law of Example 7, describing a preferred control strategy for fine-grained management of the dynamic characteristics of supercapacitor unit power output. This strategy is independent of or can be superimposed on the DICB control law, focusing on controlling the rate of power change.

[0114] In this embodiment, the step of coordinating the regulation of the energy storage unit further includes a dual-slope release law:

[0115] Step S710: Within a millisecond-level time window, the rate of change of active power of the supercapacitor unit is limited to not exceeding a preset fast power slope upper limit; and within the immediately following second-level time window, the rate of change of active power of the supercapacitor unit is limited to not exceeding a preset slow power slope upper limit.

[0116] The motivation for this step lies in the significant stratification of electromagnetic transient events in port shore power systems over time. In the initial stage of an event, such as the current spike generated at the moment of transformer closing, the duration is on the order of milliseconds, requiring the energy storage system to suppress it with extremely fast response speed and maximum power change rate to prevent a sharp drop in instantaneous voltage. However, in the middle and later stages of the event, such as the decay tail of the inrush current or the gradual change in the power feedback from the quay crane, the dynamic process slows down to the order of seconds. Based on this physical characteristic, the dual-slope release law divides the power response of the supercapacitor into two stages:

[0117] Fast response phase: Within the initial time window after the event occurs, for example, t∈[t_0,t_0+T_ms], where t_0 is the start time of the event and T_ms is the length of the millisecond-level window (e.g., 50-200ms), the rate of change of active power is allowed to reach a relatively high upper limit, i.e., |dP_supercap(t) / dt|≤s_fast.

[0118] Slow Support Phase: Within the immediately following second-level time window, for example, t∈[t_0+T_ms,t_0+T_ms+T_s], where T_s is the length of the second-level window (e.g., 1-5s), the rate of change of active power is limited to a low upper limit, i.e., |dP_supercap(t) / dt|≤s_slow. In this way, the power output of the supercapacitor can be matched to the actual demand of the disturbance, providing full support in the early stages of the event and a smooth transition in the later stages. This avoids unnecessary drastic power fluctuations, reduces stress on power electronic devices, and allows time for the intervention of other slow-regulation resources (such as lithium batteries or generators).

[0119] Step S720, wherein the upper limit of the fast power slope and the upper limit of the slow power slope are determined based on the minimum of the following physical constraints: an electrical power limit determined by the maximum current of its DC-DC converter and the minimum operating voltage of the supercapacitor; an energy release limit determined by the allowable voltage swing, total capacitance and corresponding time window length of the supercapacitor unit; and a thermal power limit derived from the thermal constraints of the device.

[0120] This step clarifies that the slope limits s_fast and s_slow are not arbitrarily set, but are determined by the device's inherent, insurmountable triple physical boundaries.

[0121] Electrical power limit: This is determined by the hardware specifications of the DC / DC converter that interfaces with the supercapacitor. Its instantaneous output power cannot exceed P_electric_limit = I_max. V_sc(t) is the maximum current that the converter semiconductor device (such as IGBT) can carry, and V_sc(t) is the real-time voltage of the supercapacitor. When calculating the most conservative upper limit of the slope, V_sc(t) is usually taken as the minimum operating voltage V_min.

[0122] Energy release limit: This is determined by the energy storage characteristics of the supercapacitor itself. Power is the rate of change of energy with respect to time. To continuously release energy at a certain power slope within a time window T_window, it inevitably leads to a reduction in the energy stored within the supercapacitor and a drop in its terminal voltage. This limit ensures that the power slope will not be so large as to drain the supercapacitor's voltage below V_min within the window period. Its approximate calculation is: P_energy_limit ≈ ΔV_allow C_sc / T_window; where: ΔV_allow is the allowable voltage drop, and C_sc is the total capacitance.

[0123] Thermal power limit: This is determined by the device's heat dissipation capacity and long-term operational reliability. Large currents flowing through the converter and the equivalent series resistance (ESR) of the supercapacitor generate Joule heat. This limit, calculated by a thermal model, is the power level at which the device can operate continuously for extended periods or at specific pulse widths without causing critical components to exceed their temperature limits; this is P_thermal_limit.

[0124] Ultimately, the preset upper limit for the fast power slope, s_fast, is determined by taking the minimum value of the above three limits on a millisecond-level time scale, that is:

[0125] s_fast=min(P_electric_limit,P_energy_limit_ms,P_thermal_limit_pulse).

[0126] Similarly, the upper limit of the slow power slope s_slow is also determined by the minimum of these three limits on a second-scale time scale, i.e., s_slow=min(P_electric_limit,P_energy_limit_s,P_thermal_limit_continuous), where the thermal constraints on a second-scale time scale are usually more stringent than those on a millisecond pulse.

[0127] In a preferred embodiment, the dual-slope release law is also prioritized in conjunction with the energy floor locking mechanism. When high-risk scenarios such as overlapping events are detected, the control system prioritizes ensuring a second-level energy floor (i.e., ensuring sufficient energy to cope with the mid-to-late stages of the event). Under this premise, it then allocates and determines the millisecond-level fast release limit based on the margin. After the risk is eliminated, the system instructs the supercapacitor to replenish energy at a slow slope s_slow until it returns to the normal steady-state voltage level.

[0128] Example 11: This example aims to illustrate in detail the practical application of the Energy Deadline Feasibility Release Criterion (EDE) in Examples 5 and 6 through a specific numerical calculation process.

[0129] In this embodiment, a specific scenario and a set of parameters are set, and the calculation and decision-making process of the EDE criterion is demonstrated step by step.

[0130] Suppose that at a certain moment, the port shore power system detects a transformer closing event, and using the method in Example 4, the impact power curve of this event is discretized into a series of cutoff-period energy requests. We select the i-th energy request for analysis.

[0131] Given the properties of the input parameter energy request i:

[0132] Energy requirement energy_i = 50 kJ. Deadline_i = 20 ms. System state and risk characteristics (both normalized): The maximum energy that the lithium battery can safely release within deadline_i is available_release_i = 600 kJ.

[0133] When the request occurs, the system's energy floor margin risk factor, floor_deficit_i, is 0.15. This value indicates that the margin between the supercapacitor's available energy and its corresponding energy floor is small, posing a certain risk.

[0134] When the request occurs, the phase misalignment risk factor phase_penalty_i = 0.10. This value indicates that the nonlinear superposition risk of the current event is at a low to medium level.

[0135] The request's slope risk factor, slope_ratio_i, is 0.80. This value indicates that the power ramp-up rate of the request is very fast, approaching the safety limit of lithium batteries, and is considered a high-risk characteristic.

[0136] Offline calibration parameters for EDE criterion:

[0137] Slope risk weight a1 = 0.5. Energy floor margin risk weight a2 = 1.2. Phase misalignment risk weight a3 = 0.8. Release threshold is 1.

[0138] The EDE criterion calculation process is as follows:

[0139] Step A: Calculate the baseline load ratio (baseline_load_ratio_i)

[0140] The assessment evaluates the proportion of the requested energy to the lithium battery's processing capacity. The calculation formula is: baseline_load_ratio_i = energy_i / available_release_i. Substituting the values:

[0141] The baseline_load_ratio_i = 50kJ / 600kJ = 0.0833 indicates that the absolute energy requirement of this request accounts for only about 8.33% of the short-term available energy of the lithium battery, which is a relatively light burden from an energy perspective.

[0142] Step B: Calculate the integrated risk coupling term (Γ_i)

[0143] This step integrates the multi-dimensional dynamic risks of the assessment request. The calculation formula is: Γ_i=a1 slope_ratio_i+a2 floor_deficit_i+a3 Substituting the value into phase_penalty_i: Γ_i = (0.5) 0.80)+(1.2 0.15)+(0.8 0.10) Component calculation: Γ_i = 0.40 (from slope risk) + 0.18 (from floor margin risk) + 0.08 (from phase misalignment risk) Summing up to: Γ_i = 0.66 This result shows that although the energy share is small, the overall risk rating of this request is 0.66, with the most significant risk contribution coming from its steep power slope (contributing 0.40).

[0144] Step C: Make a decision by applying inequality criteria.

[0145] This step integrates load and risk to make a final decision. The criterion is: baseline_load_ratio_i + Γ_i ≤ 1. Substituting the values: 0.0833 + 0.66 = 0.7433. The comparison result: 0.7433 ≤ 1, the inequality holds.

[0146] Since the calculated result of 0.7433 is not greater than the release threshold of 1, the system determines that the i-th energy request can be approved and generates an instruction to route it to the lithium battery energy storage unit for consumption.

[0147] This numerical case clearly demonstrates the decision-making logic of the EDE criterion. It doesn't simply base its judgment on the requested energy level or a single risk indicator. In this example, although the requested power slope risk is high (slope_ratio_i=0.8), its absolute energy proportion is relatively small (baseline_load_ratio_i=0.0833), and the overall system safety margin is still acceptable at the time of occurrence (floor_deficit_i and phase_penalty_i are both small). The sum of these three factors, after comprehensive evaluation, does not reach the safety boundary. Therefore, the system intelligently identifies this as a sharp but brief impact, which the lithium battery is capable of handling without jeopardizing its own safety and system stability. This avoids indiscriminately assigning all high-slope requests to the supercapacitor, thereby maximizing the utilization of the lithium battery, reducing dependence on expensive supercapacitors, and demonstrating the effect of resource synergistic optimization.

[0148] In contrast, in another scenario, the requested energy and slope remain constant, but the system's safety margin is very low when it occurs (e.g., floor_deficit_i = 0.5), then the combined risk coupling term Γ_i becomes (0.5). 0.80)+(1.2 0.5)+(0.8 0.10) = 0.40 + 0.60 + 0.08 = 1.08. At this point, the final criterion calculation result is 0.0833 + 1.08 = 1.1633, which is greater than 1, so the inequality is not true. In this case, the system will determine the request as high-risk and forcibly route it to the supercapacitor, thus prioritizing system security. This comparison further demonstrates the environmental awareness and risk adaptation capabilities of the criterion of this invention.

[0149] Example 12: This example aims to elaborate in detail the technical basis for realizing the method of the present invention—the specific construction method and calibration process of the high-fidelity electromagnetic transient model.

[0150] In this embodiment, the model construction and calibration include the following steps:

[0151] For both shore-side main transformers and ship-side receiving transformers, the models must include nonlinear magnetization branches, iron loss branches, and leakage reactance. The nonlinear magnetization branch is crucial for generating inrush current. Preferably, the Jiles-Atherton model or the Preisach hysteresis model, which accurately describes the hysteresis loop and saturation characteristics, is used. The parameters required for the governing equations of these models, such as saturation flux, coercivity, and permeability, are determined by reading the no-load test curves, short-circuit test data, and measured BH curve segments from the equipment manufacturer, combined with curves from the material handbook for fitting and completion. Furthermore, the model must include a residual magnetization initialization module to deduce the initial flux state before the current closing based on the voltage phase and demagnetization process during the transformer's previous power outage.

[0152] The upstream power grid is modeled using the Thevenin equivalent circuit, which consists of an ideal voltage source connected in series with an equivalent impedance. The equivalent impedance is set as a range parameter rather than a fixed value, depending on the different operating modes of the power grid (e.g., high / low water season, day / night load), to cover different system short-circuit capacities in subsequent offline sweeps. The transmission cables connecting the transformer and the grid connection point, depending on their length and voltage level, can be modeled using a lumped-parameter π-type equivalent circuit or a distributed-parameter line model that reflects the traveling wave process. The model parameters are determined by the cable manual and field tests.

[0153] For the Static Var Compensator (STATCOM), the quay crane feedback unit, and the interface converters for supercapacitors and lithium batteries, average models are used to accurately reproduce their external VI characteristics while ensuring computational efficiency. For example, the STATCOM model is a controlled voltage source, and its controller includes fast reactive power control and voltage support logic. The quay crane feedback unit model is an active rectifier, and its upper limit of feedback power can be used as an external control input. The supercapacitor and lithium battery are each connected to the DC bus through a bidirectional DC / DC converter, and the converter model needs to incorporate precise current limiting and power ramp-up rate constraints.

[0154] After the model is built, it needs to be calibrated to ensure it accurately reflects the dynamic behavior of the real physical system. Specifically, several typical real-world closing and inrush current waveform data can be collected from the field. Then, initial boundary conditions are set in the simulation model to be as consistent as possible with those during the actual event, and multi-objective fitting is performed on the model. The goal of the fitting is to minimize the combined error between the simulated waveform and the real waveform on the following four key indicators:

[0155] Peak current: The relative error between the first peak value of the simulated inrush current and the actual value.

[0156] Impact energy: The relative error between the simulated inrush flow and the actual energy before it decays to a certain level (e.g., 10%).

[0157] Tail-end time constant: The relative error between the time constant of the simulated inrush flow exponential decay segment and the actual value.

[0158] Harmonic content: The relative error between the simulated total harmonic distortion (THD) of the inrush flow and the actual value. This comprehensive error is minimized by adjusting parameters with high uncertainty in the model (such as core material parameters and equivalent impedance). When the relative errors of all indicators meet the preset threshold requirements (e.g., peak value and energy error not exceeding 10%, total harmonic distortion not exceeding 15%), the model can be considered calibrated, possessing sufficient fidelity, and can be used for subsequent offline domain sweep analysis.

[0159] Example 13 describes an overlapping event risk prediction method that differs from the HEAC online projection mechanism in Example 4 and can be used as an alternative or parallel option. This method is based on Markov chain theory and uses statistical modeling of historical data from port operations to predict the probability of high-risk overlapping events occurring within a short future time window.

[0160] A complete berthing and unberthing operation process can be discretized into a finite-state Markov chain. For example, this state space can be defined as eight states: S1 - Approaching berth; S2 - Mooring completed and waiting; S3 - Shore power synchronization and preparation; S4 - Shore power connection; S5 - Quay crane operation (this state can be further subdivided into sub-labels, such as active / inactive feedback); S6 - Thruster / side thruster start / stop window; S7 - Standby and transfer interval; S8 - Power outage and unberthing. Then, by analyzing historical connection and switching event records, quay crane operation logs, thruster start records, and grid connection point synchronization measurement data, a set of rules or a classifier can be established to accurately map these multi-source, continuous observation data sequences to the aforementioned discretized state space. For example, the combined observation of detecting a closing waveform and a shore power trolley position signal indicating the berth can be mapped to state S4.

[0161] After converting historical data into a state sequence, a state transition count matrix is ​​obtained by counting the frequency of state pairs occurring at adjacent time points. For example, the number of transitions from state S4 to state S5 is counted. Based on this count matrix, the maximum likelihood estimation method is used to calculate the state transition probability matrix P of the discrete-time Markov chain, where the element P_ij represents the probability of transitioning from state i to state j. To handle the data sparsity problem, Laplace smoothing or Good-Turing smoothing is preferably used when calculating the probabilities. In some alternative implementations, if the dwell time distribution of some states is observed to be significantly non-geometric (i.e., memoryless), the model can be upgraded from a discrete-time Markov chain (DTMC) to a semi-Markov model, where state transitions depend not only on the transition probabilities but also on the dwell time distribution in that state, which can be represented by an empirical distribution or a phase-type distribution fitted from historical data.

[0162] To improve the prediction accuracy of the model, the Markov chain can be constructed as non-homogeneous, meaning its state transition probabilities are not fixed but depend on external contextual variables. For example, the operating cycle time and state transition probabilities of large container ships and small bulk carriers are obviously different; daytime operating efficiency may also differ from nighttime efficiency. Therefore, the transition probabilities can be made to explicitly depend on contextual variables such as time of day (day / night), berth, ship type, and tidal window. This can be achieved by binning the data according to context and estimating a transition matrix separately for each bin; or by using a generalized linear model (GLM) with contextual variables as independent variables in a regression model to predict the transition probabilities.

[0163] After establishing the above model, online risk prediction can be performed. For example, when the system enters state S4 (shore power connection), the state transition matrix can be used to calculate the joint probability that the system will simultaneously enter the feedback active sub-state of state S5 and the thruster start / stop window of state S6 within the next N steps (e.g., the next 5 minutes). This calculated probability is the predicted probability value of the short-term overlapping event. When this predicted probability exceeds a preset risk threshold, the system will trigger an early warning and can set a more conservative energy floor target accordingly, or pre-charge the energy storage system in advance. This probability prediction result can also be directly used in opportunity-constrained programming problems. For example, when determining the energy floor target E_floor, the constraint P{overlapping energy demand ≤ E_floor} ≥ 1 - ε is applied, where ε is the acceptable risk level.

[0164] Example 14: This example aims to elaborate on a key prerequisite step in the robustness analysis and security certificate generation of this invention: the method for constructing the parameter uncertainty set. It is understood that a precise and reasonable set describing all sources of uncertainty in the system is the cornerstone for all subsequent offline simulation analyses and online robust control to arrive at reliable conclusions. This example provides specific construction basis and procedures for the parameter uncertainty set mentioned in Examples 2, 9, and others.

[0165] In one specific embodiment of the present invention, the process of constructing the parameter uncertainty set includes the following steps:

[0166] First, it is necessary to systematically identify the key parameters that affect the electromagnetic transient processes of the port shore power system and whose values ​​fluctuate or cannot be accurately determined. In this embodiment, these parameters are identified and classified as follows:

[0167] Event-driven uncertainties are mainly related to the triggering conditions of impact events. For example, the phase angle at which the three-phase contacts close at the moment of transformer closing is not a precisely controllable constant, but rather exhibits millisecond-level jitter and asymmetry, i.e., closing angle jitter. In addition, the residual magnetization level and polarity left by the transformer after the last power outage also have strong randomness.

[0168] The inherent uncertainties of equipment mainly stem from the manufacturing tolerances of electrical equipment. For example, the actual resistance of a pre-inserted resistor used to suppress inrush current will fluctuate within a specific error range (e.g., ±5%) around its nominal value.

[0169] The uncertainty of system operation status is mainly related to the operation mode of the upstream power grid. For example, the equivalent impedance or short-circuit capacity of the upstream power grid will fluctuate within a large range as generator units in the grid start and stop, load levels change, and network topology changes.

[0170] After identifying the sources of uncertainty, it is necessary to quantify and model each type of parameter mathematically.

[0171] Quantification of random parameters: For parameters with strong randomness, such as closing angle jitter and residual magnetization, the quantification process is as follows: First, by analyzing a large amount of historical event data collected from a high-precision field recording device or a synchronous phasor measurement unit (PMU), the actual occurrence sequence of these parameters is extracted. Then, statistical analysis is performed on this sequence to calculate its mean, variance, skewness, and other statistical moments. Based on this, preferably, a probability density function that best describes the distribution characteristics of the parameter is fitted. For example, closing angle jitter can usually be approximated as a Gaussian distribution with the expected closing angle as the mean. Finally, the uncertainty of the parameter is modeled as a specific probability distribution. In some optional implementations, if the historical data is sparse or the distribution characteristics are not obvious, non-parametric methods such as kernel density estimation can also be used for modeling, or a bounded interval covering the vast majority of historical data (e.g., the mean plus or minus three standard deviations) can be used to conservatively describe its uncertainty.

[0172] For the quantification of interval parameters: For parameters such as the short-circuit capacity of the upstream power grid and equipment manufacturing tolerances that vary within a certain range over a specific time scale, the quantification process is as follows: by consulting system analysis reports under different operating modes provided by the power grid planning department, or product data sheets provided by equipment suppliers, the upper and lower limits of their variation are directly determined. For example, according to the power grid plan, the short-circuit capacity of a substation connected to a port can reach up to 1500 MVA under high load conditions during the high-water season, and may drop to as low as 500 MVA under low load conditions during the dry season. The uncertainty of this parameter is then modeled as a closed interval [500 MVA, 1500 MVA].

[0173] All quantified and modeled uncertainty parameters are combined into a formalized mathematical set, namely the parameter uncertainty set, usually denoted as Θ. This set is a multi-dimensional parameter space, where each dimension corresponds to an uncertainty parameter. For example, a simplified parameter uncertainty set can be defined as: Θ={θ|θ=(θ_angle,θ_resistor,θ_scc,θ_flux) and θ_angle∈N(μ_a,σ_a^2), θ_resistor∈[R_nom (1-tol),R_nom (1+tol)],θ_scc∈[SCC_min,SCC_max],θ_flux∈[Flux_min,Flux_max]}. This formally defined set will serve as the input for subsequent offline domain scanning simulations (such as in Example 2). The simulation program will sample within this set to ensure that the generated safety certificate covers all possible combinations of operating conditions.

[0174] Example 15: This example illustrates the complete calculation process from the HEAC segment (power upper bound) to the quantile energy floor baseline (statistical energy upper bound).

[0175] According to one aspect of this application, a method for multi-timescale capacity configuration of a port energy storage system further includes:

[0176] Step S1010: Perform time integration on the HEAC segment in the security certificate library to obtain the cumulative demand curve characterizing the cumulative energy demand under the worst-case scenario.

[0177] Specifically, HEAC fragment A mix (t) describes the worst-case power upper bound at the boundary of a specific set of parameter uncertainties (or a specific operating condition). This function is typically piecewise affine or piecewise linear. The time integral, which is the definite integral of this function from the event's initial time 0 to a future time t, is calculated as: E cum (t)= [0,t]A mix (τ)dτ. Since A mix (t) is a piecewise function, and its integral can be calculated analytically piecewise, resulting in higher computational efficiency. The resulting cumulative demand curve E cum The physical meaning of (t) is the total energy required up to time t to completely absorb the impact under the corresponding specific working condition.

[0178] Step S1020: Sampling of the cumulative demand curve based on the set of parameter uncertainties.

[0179] E cum (t) represents only the result for one operating condition. To obtain a statistical upper bound that covers all uncertainties, scenario sampling is required. This step involves generating a large number of parameter combination scenarios from the parameter uncertainty set Θ constructed in Example Fourteen. Specifically, this sampling process can be carried out in one of the following ways:

[0180] Latin hypercube sampling (LHS): This method divides the marginal distribution of each parameter into multiple equally probable intervals and ensures that only one sample is drawn from each interval. Compared to simple random sampling, LHS can achieve a more uniform and comprehensive coverage of the entire multidimensional parameter space with fewer samples (e.g., hundreds to thousands of samples), thereby significantly reducing the computation time for subsequent electromagnetic transient simulations while maintaining accuracy.

[0181] Monte Carlo random sampling (MC): This method generates a parameter combination vector by independently sampling from the probability distribution (e.g., Gaussian distribution) or within its interval (e.g., uniform distribution) of each parameter. This method is simple and intuitive, but to achieve sufficient coverage, it typically requires a very large number of samples (e.g., tens of thousands of samples). For each sampled parameter scenario j, the steps in Example 2 are repeated to calculate its specific HEAC fragment. \mathcal{A}_{mix,j}(t) And by integrating it, we obtain its own cumulative demand curve E^cum_j(t).

[0182] Step S1030: Take quantiles according to the preset confidence level to generate the quantile energy floor baseline curve.

[0183] After obtaining a family of cumulative demand curves {E^cum_1(t), E^cum_2(t), ..., E^cum_N(t)} containing N scenarios through scenario sampling, this step aims to extract a statistically significant upper boundary from this family of curves. The calculation process is as follows: at each discrete time point t_k, collect the cumulative demand values ​​of all N scenarios at that time point, forming a numerical set {E^cum_1(t_k), E^cum_2(t_k), ..., E^cum_N(t_k)}. Then, sort this set of N values ​​in ascending order. Quantile calculation, that is, according to a preset confidence level q (for example, q = 0.95, representing a 95% confidence level), extract the q-th quantile from the sorted set. The values ​​at N locations (linear interpolation performed if necessary) are used to determine the quantile energy floor baseline value E_floor_base(t_k) at time t_k. This process is repeated for all time points t_k of interest, and the resulting series of E_floor_base(t_k) values ​​are concatenated to form the final quantile energy floor baseline curve. This curve statistically guarantees a 95% probability that the worst-case cumulative energy demand will not exceed it.

[0184] Understandably, the choice of the pre-set confidence level q represents a trade-off between security and economy. Choosing a higher confidence level (e.g., q=0.99 or 0.999) will result in a higher energy floor baseline, making the system more secure and resilient to more extreme and rarer combinations of parameter events, but correspondingly requiring a larger (and therefore more expensive) energy storage system. Conversely, choosing a lower confidence level can reduce the required energy storage capacity and investment costs, but requires accepting a known, higher probability of failure. The choice of this value is typically determined by the port's operational safety level, grid specifications, and the project's return on investment analysis.

[0185] This invention fundamentally solves the technical problem of accurately quantifying and assessing the superposition risk of asynchronous concurrent transient impacts in the background technology by constructing and applying the HEAC safety certificate library, achieving a precise and tight definition of the system safety boundary. This effect is achieved by introducing an innovative hybrid electromagnetic arrival calculus theory adapted to electromagnetic transient analysis. Specifically, unlike traditional methods that oversimplify or conservatively superimpose transient events, this invention utilizes the phase penalty function as a core mathematical tool to quantitatively characterize for the first time the strong nonlinear and non-commutative superposition cost caused by microsecond to millisecond-level phase misalignment between events such as transformer inrush current and quaybridge energy feedback. Based on this, by finding the minimum upper envelope over all possible phases, a HEAC power upper bound that can cover the worst case while avoiding over-conservatism is generated. Furthermore, by sampling the parameter uncertainty set for scenarios and taking quantiles, this power upper bound is transformed into a statistically significant quantile energy floor baseline. For the specific scenario of shore power in ports, this precise boundary definition avoids serious over-allocation of energy storage capacity due to simple arithmetic superposition, saving millions or even tens of millions of yuan in initial investment costs. On the other hand, it overcomes the shortcomings of traditional power flow simulation in failing to capture transient risks, providing the system with a scientific and reliable safety margin, eliminating the risk of voltage collapse and shutdown under extreme operating conditions due to capacity under-allocation, and achieving a high degree of unity between safety and economy.

[0186] This invention addresses the shortcomings of background technologies in energy storage control, namely the lack of foresight and the inability to provide deterministic safety guarantees, by applying a demand integral barrier-guaranteed control law. This provides supercapacitor units with a rigid safety guarantee that their energy floor cannot be breached. This effect is achieved by innovatively transforming the control domain from the traditional state trajectory domain to the energy integral domain. Specifically, the control law constructs a barrier function B(t) directly related to the current available energy of the supercapacitor and the predicted cumulative energy demand in the future. Its core control objective is to maintain the non-negativity of B(t). By solving a control solution that maintains the barrier function in accordance with a specific preserving inequality, DICB can achieve pre-emptive defense against future energy shocks, rather than the delayed response of traditional control strategies after adverse consequences such as voltage drops. In scenarios like port power grids where shocks are unpredictable and have severe consequences, this forward-looking control mechanism with rigorous mathematical guarantees is crucial. It ensures that the supercapacitor, which serves as the last line of defense for system safety, always has sufficient energy reserved for its core, millisecond-level rapid response tasks, greatly improving the power supply reliability of critical port operations (such as synchronous hoisting of quay cranes) and avoiding controller downtime or operation interruption caused by instantaneous power shortages.

[0187] This invention addresses the technical problems of coarse-grained hybrid energy storage control strategies and the inability to differentiate power requests with varying risk characteristics in the prior art by designing and applying the EDE release criterion. This achieves refined protection of lithium batteries and optimization of overall system operating costs. This effect stems from its innovative multi-dimensional risk perception and routing decision-making mechanism. Unlike traditional blind division of labor methods based on high- and low-frequency filtering, the EDE criterion comprehensively and quantitatively evaluates each discrete energy request from three key dimensions before making a scheduling decision: power slope (impact on battery life), energy floor margin (threat to current system safety), and phase misalignment (inherent nonlinear risk of the event). Only requests deemed low-risk by this combined criterion are released to the lithium battery. For port hybrid energy storage systems, this means that seemingly mild but potentially dangerous requests occurring during vulnerable moments or containing severe nonlinear risks can be accurately identified and processed by supercapacitors, effectively avoiding potential damage to lithium batteries and extending their cycle life by more than 20%, according to calculations. This intelligent energy routing maximizes the role of lithium batteries as a low-cost energy source while ensuring safety, significantly reducing the levelized cost of electricity (LCOE) over the entire system's lifecycle.

[0188] This invention addresses the common problem in existing energy storage configurations—the disconnect between static, offline models and dynamic, changing actual operating conditions—by introducing a closed-loop tuning and updating mechanism for the safety certificate repository. This endows the system with self-learning and adaptive capabilities. This effect is achieved through a feedback loop of monitoring-comparison-triggering-updating. During actual operation, the system continuously collects key verification indicators and compares them with the theoretical boundaries of the safety certificate repository. Once it detects that real-world behavior (such as an unexpected surge) crosses the theoretical boundary, or that the theoretical boundary is too conservative, the system automatically triggers the certificate update process. It traces the parameters that caused the deviation (such as an unexpectedly low grid impedance), adjusts the parameter uncertainty set accordingly, and then re-executes the entire offline simulation and certificate generation process, resulting in a more accurate and realistic new certificate repository. This adaptive capability is particularly important for the constantly evolving scenario of ports (e.g., adding berths, replacing transformers, and upgrading the upstream power grid). It ensures that the safety and economy of the energy storage system will not decline with changes in the external environment, guarantees the long-term effectiveness and advancement of the solution over a lifespan of ten years or even longer, avoids costly manual remodeling and tuning, and achieves true intelligent operation and maintenance.

[0189] This invention discloses a multi-timescale capacity configuration method for port energy storage systems, aiming to solve the problems of inaccurate quantification of transient impact superposition risks and coarse-grained collaborative control of hybrid energy storage in existing technologies. The method includes: an offline stage, based on an innovative Hybrid Electromagnetic Arrival Algorithm (HEAC) theory, using electromagnetic transient simulation and upper envelope extraction, and applying a phase penalty function to define the asynchronous superposition cost, constructing a safety certificate library containing HEAC segments; an online stage, real-time matching of early features of synchronous measurement data from grid connection points with the safety certificate library, generating a time-based energy floor curve for real-time control; finally, applying a Demand Integral Barrier Guarantee (DICB) control law to provide deterministic safety guarantees for supercapacitors, and simultaneously applying the Energy-Deadline Feasibility Release (EDE) criterion for risk perception energy routing of lithium batteries, collaboratively controlling energy storage units based on the energy floor curve, and finally determining the operational control scheduling table and capacity configuration result list. This invention achieves safe, economical, and refined configuration and scheduling of port energy storage systems across multiple timescales.

[0190] The Hybrid Electromagnetic Arrival Algorithm (HEAC) theory proposed in this invention is innovative in that it solves the fundamental problem that traditional arrival curve theory in network calculus cannot be directly applied to the electromagnetic transient analysis of power systems. Traditional arrival curve theory is based on the assumption of first-in-first-out (FIFO) data packet networks and a fixed service rate. However, the superposition of electromagnetic transient events such as transformer inrush current and quay bridge energy feedback is a highly nonlinear physical process strongly correlated with the phase angle of the event occurrence, and does not satisfy simple linear superposition or subadditivity. This invention creatively introduces a phase penalty function... (Δ) represents the first mathematical quantification of the nonlinear superposition cost caused by phase misalignment. Its physical meaning can be understood as an additional power / energy compensation necessary to ensure that the simple arithmetic of two asynchronous events can safely encompass their true, complex physical superposition result. This enables HEAC theory to generate more compact and formally bounded safety certificates than traditional transient analysis methods (such as time-domain simulation superposition), marking a crucial step from stochastic, discrete simulation analysis to deterministic, boundary-oriented theoretical analysis.

[0191] The Energy-Deadline Feasibility Release Criterion (EDE) proposed in this invention achieves multi-dimensional integrated innovation in energy storage resource scheduling. The coupling of three dimensions—slope, energy floor margin, and phase misalignment—is necessary because single-dimensional decision-making has significant drawbacks: focusing solely on energy magnitude ignores the severity of power surges; focusing solely on power slope fails to assess the system's overall safety reserves at that time. This invention recognizes that whether a lithium battery can safely respond to a request is the result of the combined effects of its own dynamic capabilities (slope), the system's current safety state (floor margin), and the inherent risk of the event (phase). Using a linear combination is an efficient and interpretable approximation of this multi-dimensional decision space, possessing engineering sufficiency while meeting the computational speed requirements of real-time control. It differs fundamentally from the earliest deadline-first (EDF) algorithm in traditional real-time scheduling: EDF is purely time-driven, only concerned with whether the task can be completed before the deadline, without considering the cost of task execution; EDE, on the other hand, is risk-driven. While meeting time constraints, its core is to assess the safety cost of task execution, thereby making decisions to protect high-value assets (lithium batteries) and ensure system safety.

[0192] The Demand Integral Barrier Guarantee Control Law (DICB) proposed in this invention represents a fundamental innovation compared to the traditional Control Barrier Function (CBF) method. Traditional CBF typically defines a safety set in the system's state space (such as the instantaneous value trajectories of voltage and current), while the DICB of this invention achieves a transformation from the state trajectory domain to the energy integral domain. That is, the barrier function B(t) is directly related to the available energy of the supercapacitor. E_{supercap}(t) The future cumulative energy demand of the event \widehat{E}^{cum}(t) The advantages of this transformation are: 1. It directly manages the root cause of system instability—energy imbalance—rather than merely its manifestations (such as voltage drops), making its control effect more fundamental. 2. By pre-integrating, it uses the total amount of future shocks as the basis for current decisions, providing foresight and enabling earlier responses. 3. The mathematical rigor of its guaranteed control lies in the fact that as long as the maintenance inequality holds, the risk of the energy storage unit's usable energy being depleted by future shocks can be theoretically eliminated, providing a deterministic mathematical guarantee for system safety.

[0193] According to one aspect of this application, the dual-slope release law described in Embodiment 8 may further include a priority and dominance rule: when the system identifies a high-risk scenario of overlapping events (e.g., shore power closing and shore bridge full-load feedback occurring simultaneously) through access and switching event records or short-term overlap probability prediction (as described in Embodiment 13), the dominance rule is activated. The specific implementation of this rule is as follows: the control system will prioritize ensuring energy balance within a second-level time window, i.e., the constraint level of the DICB control law described in Embodiment 7 is raised to the highest level to ensure the system has sufficient energy to weather the entire event. Under this premise, the upper limit of the fast power slope for the millisecond-level time window is... s_{fast} It will no longer be a fixed value, but a variable that is dynamically adjusted based on the energy margin calculated by DICB. The smaller the margin, the lower the energy margin. The lower the upper limit of s_{fast} is, the more likely the risk event will be (e.g., when the inrush current is detected to have decayed to a steady-state value of 5s_{slow}). To constrain the ramp rate, the supercapacitor is instructed to absorb power from the grid smoothly until its terminal voltage recovers to the steady-state operating point.

[0194] According to one aspect of this application, the objective function and constraint set of the joint optimization problem described in Embodiment Nine can be specifically constructed as follows:

[0195] Mathematical expression of the objective function: The joint optimization objective function aims to minimize the total lifecycle cost, and its specific form is as follows:

[0196] min(w_cap Cost_cap(P_sc,E_sc,P_li,E_li)+Σ_[t∈T](w_age Cost_age(I_li(t),ΔSOC_li(t))+w_pen Cost_penalty(V(t),P_rev(t)))).

[0197] Where: P_sc, E_sc, P_li, and E_li are the power and energy capacity of the supercapacitor and lithium battery to be optimized, respectively; Cost_cap(·) is the annualized investment cost function related to capacity; T is the total optimization time; Cost_age(·) is the aging cost function of the lithium battery, which is a function of its real-time current I_li(t) and cycle depth ΔSOC_li(t); Cost_penalty(·) is the power quality violation penalty function, which is the penalty for exceeding the soft constraint of the grid connection point voltage V(t) and the reverse power flow power P_rev(t); w_cap, w_age, and w_pen are the weighting coefficients of each cost.

[0198] The complete definition of the constraint set: The constraint set ST of this optimization problem includes at least the following:

[0199] Power balance constraint: ΣP_output(t) + P_load(t) = P_grid(t)

[0200] HEAC certificate upper bound constraint: ΣP_output_max(t)≥Apredicted(t);

[0201] DICB preserves the inequality constraint: B[k+1]-(1-k) T_s) B[k]≥0;

[0202] EDE release criterion constraint: (1-z_i) energy_i=0,ifbaseline_load_ratio_i+Γ_i>1, where z_i is a binary variable indicating whether request i should be allocated to the lithium battery.

[0203] Physical constraints for each device: P_min≤P(t)≤P_max, SOC_min≤SOC(t)≤SOC_max, etc.

[0204] Solution Algorithm and Compilation Process: This problem is essentially a mixed-integer nonlinear programming (MINLP) problem. For long-term capacity configuration, heuristic algorithms such as genetic algorithms or particle swarm optimization can be used. For short-term runtime scheduling, it can be relaxed or decomposed into a series of quadratic programming (QP) problems solvable in minutes. The compilation process refers to converting the decision variables output by the optimizer (e.g., P_sc = 1.2MW at t=5s) into specific instructions that the underlying hardware can execute. For example, the floating-point number 1.2MW is converted into a 16-bit integer (e.g., 32768) through table lookup or formula conversion, and then written to a specific register address of the PLC or DSP controlling the supercapacitor PCS (power conversion system) via communication protocols such as Modbus-TCP or IEC61850, thereby completing the control closed loop.

[0205] According to one aspect of this application, the verification and update process described in Embodiment 10 can be specifically described as follows: the scenario playback method can preferably use a hardware-in-the-loop (HIL) or software-in-the-loop (SIL) test platform to inject the collected historical disturbance data as an excitation signal into the real-time simulation model of the controlled energy storage system.

[0206] The specific calculation of the verification index set, for example, the grid connection point voltage deviation index can be specifically calculated as the integral absolute value of the deviation between the effective voltage value and the rated voltage over the entire event window. The parameter uncertainty set adjustment method is as follows: when a boundary crossing caused by excessively low upstream grid impedance is detected, the system will automatically identify the operating condition label corresponding to the event and locate the defined interval [SCC_min, SCC_max] of the parameter uncertainty set Θ regarding the grid equivalent impedance θ_scc. Subsequently, the lower limit SCC_min of this interval is reduced by a preset step size (e.g., 5% of the original lower limit) to form an updated uncertainty set, which is then used as input to trigger the regeneration of the certificate.

[0207] The preferred embodiments of the present invention have been described in detail above. However, the present invention is not limited to the specific details in the above embodiments. Within the scope of the technical concept of the present invention, various equivalent transformations can be made to the technical solutions of the present invention, and these equivalent transformations all fall within the protection scope of the present invention.

Claims

1. A method for multi-time scale capacity configuration of a port energy storage system, characterized in that, The method comprises the following steps: The off-line establishment of the security certificate library comprises the following steps: The online matching and grid point synchronous measurement data and the security certificate library are used to generate a time energy floor curve for real-time control; The application of the demand integral obstacle guarantee type control law and the energy cut-off period feasibility release criterion cooperatively regulates the energy storage unit according to the time energy floor curve, determines an operation period control schedule and a capacity configuration result list; The off-line establishment of the security certificate library comprises the following steps: The phase penalty function is calibrated through the following steps:

2. The port energy storage system multi-time scale capacity configuration method according to claim 1, characterized in that, The application of the demand integral obstacle guarantee type control law cooperatively regulates the energy storage unit, comprising: The online matching and grid point synchronous measurement data and the security certificate library are also used to generate a phase-corrected cumulative demand curve; The demand integral obstacle function is constructed according to the available energy of the super capacitor unit, the phase-corrected cumulative demand curve and a preset error upper limit margin; The control solution for maintaining the demand integral obstacle function to satisfy the preset maintenance inequality is solved, thereby generating a power reference trajectory for the super capacitor unit to guarantee that the energy floor is not crossed.

3. The port energy storage system multi-time scale capacity configuration method of claim 1, wherein, The application of the energy cut-off period feasibility release criterion cooperatively regulates the energy storage unit, comprising: The online matching and grid point synchronous measurement data and the security certificate library are also used to discretize the impact event into a cut-off period energy request set carrying energy demand and cut-off period information; For any energy request in the cut-off period energy request set, the baseline load ratio is quantified, and the coupling factor representing the risk of the energy request is calculated from three dimensions of slope, energy floor margin and phase offset; According to the combination criterion of the baseline load ratio and the coupling factor, the energy requests that can be safely accepted by the lithium battery energy storage unit are identified, and the remaining energy requests are allocated to the super capacitor unit to form an energy routing strategy.

4. The port energy storage system multi-time scale capacity configuration method according to claim 3, characterized in that, The weight coefficient determination process of the coupling factor in the combination criterion is as follows: An EDE calibration sample set is constructed, in which, by replaying the scenario that any energy request is taken by the lithium battery energy storage unit, if the hard constraint violation is triggered, the energy request is labeled as unsafe, otherwise, it is labeled as safe; A linear programming problem is established and solved, in which, the hard constraint is that all energy requests carrying unsafe labels must be blocked, and the optimization objective is to minimize the number of energy requests carrying safe labels that are rejected, so as to obtain the weight coefficient.

5. The port energy storage system multi-time scale capacity configuration method of claim 1, wherein, Further comprising: In the playback of critical scenarios or actual operation, the calibration index set under the execution of the operation period control schedule table is collected; The calibration index set is compared with the boundary guarantee of the security certificate library, and when the boundary crossing or the error approximation error is over the limit, the certificate update is triggered; In response to the certificate update, the parameter uncertainty set related to the boundary crossing or the error over the limit is adjusted, and the step of establishing the security certificate library offline is re-executed to realize the closed-loop tuning of the security certificate library.

6. The port energy storage system multi-time scale capacity configuration method of claim 1, wherein, The step of online matching and grid synchronization measurement data with the security certificate library includes: From the millisecond level starting segment of the grid synchronization measurement data, the online early feature sequence containing high frequency voltage distortion, initial magnetization response and closing asymmetry index is extracted; Based on the online early feature sequence, set consistency matching is performed in the HEAC segment of the security certificate library to identify the most consistent HEAC segment and estimate the time phase offset; The most consistent HEAC segment and the time phase offset are combined to synthesize the time type energy floor curve.

7. The port energy storage system multi-time scale capacity configuration method of claim 3, wherein, The combination criterion specifically represents an inequality criterion: The coupling term factor of the three dimensions of slope, energy floor margin and phase misplacement is weighted and summed to form a comprehensive risk coupling term; When the baseline load ratio and the sum of the comprehensive risk coupling term are not greater than the preset release threshold, it is determined that the energy request can be safely taken by the lithium battery energy storage unit.

8. The port energy storage system multi-time scale capacity configuration method of claim 1, wherein, The steps of determining the operation period control schedule table and the capacity configuration result list include: An integrated optimization objective function is constructed, which combines the device capacity cost, lithium battery aging cost and voltage and reverse power flow violation penalty; And under the joint constraints of the certified upper bound provided by the security certificate library, the maintenance inequality of the demand integral barrier type control law and the energy-cut-off period feasibility release criterion, the optimization result set is solved; The optimization result set obtained by solving is compiled into the operation period control schedule table, and the capacity configuration result list containing the energy storage unit capacity and operation constraint parameters is derived.

Images