Modelling of Supercapacitor Banks for Power System Dynamics Studies
The paper presents accurate and simple dynamic model of a supercapacitor bank system for power system dynamics studies. The proposed model is derived from a detailed RC circuit representation. Furthermore, a complete control system of the supercapaci…
Authors: Matej Krpan, Igor Kuzle, Ana Radovanovic
1 Modelling of Supercapacitor Banks for Po wer System Dynamics Studies Matej Krpan, Student Member , IEEE, Igor Kuzle, Senior Member , IEEE, Ana Radov anovi ´ c, Student Member , IEEE, Jovica V . Milanovi ´ c, F ellow , IEEE Abstract —The paper presents accurate and simple dynamic model of a supercapacitor bank system for power system dynam- ics studies. The pr oposed model is derived from a detailed RC circuit repr esentation. Furthermore, a complete control system of the supercapacitor bank is also pr esented. The proposed model is easy to integrate in any power system simulation software and consists of only up to four standard datasheet parameters. The performance of the proposed model in grid frequency control and low-voltage ride through is illustrated on IEEE 14-bus test system in DIgSILENT PowerF actory . It is sho wn that in case of transient stability simulations the ideal (simplified) model of the supercapacitor can be used while in case of frequency contr ol the ideal representation may not always be appropriate. Index T erms —po wer system dynamics, power system simula- tion, power system modeling, power generation control I . I N T RO D U C T I O N T HE trend of increasing in verter-interfaced generation (IIG) in power systems throughout the world and subse- quent reduction of synchronous inertia has motiv ated many re- search efforts to understanding stability of low-inertia systems as well as de veloping new algorithms which enable the IIG participation in system frequenc y control and other ancillary services. Supercapacitor (SC) energy storage system (ESS) can be used both as a standalone ESS for grid support or as a com- bination with other storage systems or IIG as part of a hybrid energy (storage) system. Its high power density , in particular , as well as hundreds of thousands of charging/dischar ging cycles and fast discharge naturally make it most applicable during grid frequenc y excursions when a fast injection or absorption of acti ve po wer is necessary . Similarly , it can be used to quickly stabilize intermittent output of solar and wind generation. There are sev eral reasons for using SC systems for fast injection of high power instead of other storage devices, e.g. batteries or flywheels [1], [2]: i) SC bank can be fully charged or dischar ged in the time scale of se veral tens of seconds or faster while the rated po wer can be reached within a few milliseconds; ii) SCs have bigger power density than batteries and flywheels, i.e. a SC system of the same power rating will be much smaller than equiv alent battery or flywheel M. Krpan and I. Kuzle are with the Department of Energy and Power Systems, Faculty of Electrical Engineering and Computing, Uni versity of Zagreb, 10000 Zagreb, Croatia, e-mail: matej.krpan@fer .hr; igor .kuzle@fer .hr A. Radovanovi ´ c and J. V . Milanovi ´ c are with the School of Electrical and Electronic Engineering, University of Manchester , Manchester , UK, e-mail: ana.radov anovic@manchester .ac.uk; milanovic@manchester .ac.uk This work was supported by the EU Horizon 2020 project “CROSSBO W”, grant agreement No. 773430. system; iii) SCs can withstand significantly more (hundreds of thousands) charging/dischar ging cycles; iv) SCs hav e smaller operation and maintenance costs than batteries and flywheels. A. Literatur e survey SC technology has often been used for electric vehicle appli- cations in the past, e.g., [3], with the focus on numerical mod- elling and/or energy management. SCs are often paired with wind and/or solar generation systems for power smoothing, virtual inertial response or lo w-voltage ride through (L VR T) [4]–[7] and in these papers the SC is usually used in the DC link of voltage source con verters. On larger scale applications, SC is often used as a part of a hybrid ESS in microgrids, e.g. [8], or isolated power systems [9]–[14] for lev eling out intermittent rene wables or for grid ancillary services such as frequency or voltage support. All of the surveyed papers hav e one or more similarities: i) SC is modelled as an ideal capacitor which is not alw ays appropriate as the capacitance, and therefore the stored energy as well, varies with the applied voltage [15]; ii) SC energy storage system model is not applicable for power system dynamic studies as the supercapacitor model is either giv en in its RC/RLC form (not a block diagram with defined inputs and outputs), e.g. [12], [16] or the complete control system is only gi ven functionally on a higher level with actual details on subsystem implementation missing which makes it difficult to integrate in power system simulation software, e.g. [13]. One outlier from the revie wed literature is a simplified model by Egido et al. [14] which neglects any capacitor dynamics and is described only with an initial state-of-charge (SoC) estimation and a simplified control system which seems to agree fairly well with field measurements. Howe ver , the used disturbances were smaller than the size of a fully charged supercapacitor and the time scale was not long enough to observe the differences in time-to-discharge (i.e. when the stored energy is depleted). As will be shown in this paper, both initial SoC and the size of a disturbance impact the performance of different supercapacitor models. T o the best of our knowledge, there were no detailed studies on the adequate complexity le vel of a SC model for po wer system dynamics studies, no rigorous deriv ation of the SC bank model and no complete representation of a SC bank system with associated control that can be easily integrated in commercially av ailable power system simulation software. Hence, this paper is trying to bridge the identified gap in modeling. 2 B. Contribution A complete supercapacitor/ultracapacitor bank model should include: accurate dynamics of a SC cell, SC DC current calculation, charge/discharge control, active power and voltage/reacti ve power in verter control as well as frequency control loop. Therefore, the contributions of this paper are as follows: i) detailed analysis and comparison of different SC models with varying le vels of detail; ii) dev elopment of the appropriate models of different components of SC based energy storage system for po wer system dynamics applications (voltage and frequency control, transient stability studies); iii) dev elopment of accurate dynamic model of the SC bank with all the necessary controls that is easy to integrate in any commercial power system simulation software. I I . D E R I V A T I O N O F S U P E R C A P AC I T O R B A N K M O D E L A. Supercapacitor theory Core of the SC bank model is the SC cell. Overvie w of different SC models can be found in [17] while the state-of- the-art SC models are av ailable in [15], [17], [18]. Basically , all these models are based on RC circuit identification using impedance spectroscopy . They can be described with the same type of RC circuit consisting of three parallel sections as shown in Fig. 1. The first branch { M1 } models fast dynamics, parallel branches { M2 } model slo wer recombination phenom- ena after a fast charge or discharge and the third branch { M3 } models the long-term self-discharge phenomena [15]. The time constants of the RC circuits in the parallel branches (ranging from a minute or sev eral minutes up to an arbitrary amount of time as reported in [15], [17]) are much longer than the timescale in vestigated in this paper (up to around 30 s), hence the y can be neglected in SC models dev eloped here as will be shown in section II-B. The important characteristics of SCs for model dev elopment can be summarized as follows: i) majority of the ultracapac- itor capacitance comes from C sc ; ii) series combination of parallel branches R s 1 C s 1 – R s n C s n is actually an infinite series of these parallel groups. Howe ver , it has been shown in [15] that five groups are sufficient to obtain an accurate model; iii) capacitance C sc as well as infinite sum elements R s k , C s k are dependent on the ultracapacitor voltage u C ( t ) . This dependence is nonlinear hence the model has the time-v arying parameters. That is why an ideal capacitor representation used in many papers in the past is not always appropriate; iv) the number of parallel branches in the { M2 } group is also theoretically infinite, though it has been sho wn that between two and four branches are sufficient to achiev e accurate results [15]. In order to simplify the model further the following as- sumptions are made: i) R s is the equiv alent series, voltage dependent resistance (ESR) determined at very high frequency [15]. Considering that R s is small ( < 10 m Ω ) and its impact on the model performance is insignificant, it is typically considered to be constant; ii) temperature dependence of the parameters is neglected, i.e., the temperature is considered to be constant. The assumption is that the cooling of the system is adequate and that the system operates at room temperature. { M1 } :First branch with n parallel groups { M2 } : m parallel branches { M3 } :Leakage curren t branch i sc R s + − u R s C sc ( u C ) + − u C C s 1 ( u C ) + − u C s 1 R s 1 ( u C ) C s n ( u C ) + − u C s n R s n ( u C ) + − u sc ( t ) R p 1 C p 1 R p m C p m R leak Fig. 1. Detailed RC circuit of a supercapacitor cell Although this effect can be included in the model, it was considered that for the initial deri vation of the model for po wer system dynamic studies the temperature can be considered to be constant. Parameters of the first branch are calculated according to (1)–(3) [15]. C sc ( u C ) = C 0 + k v u C ( t ) (1) C s k = 1 2 C sc , k ∈ { 1 ...n } (2) R s k = 2 τ ( u C ) k 2 π 2 C sc (3) C 0 is the ultracapacitor capacitance at 0 V and k v is a constant expressed in F/V . τ ( u C ) is another experimentally determined parameter (it has a dimension of time) that can be expressed as function changing linearly with the voltage u C : τ ( u C ) = τ 0 + k τ u C ( t ) [15]. Howe ver , it can also be approximated by (4) [15]: τ ( u C ) ≈ 3 C sc ( R dc − R s ) , (4) where R dc is the resistance experimentally obtained at v ery low frequencies (essentially DC). Naturally , R dc > R s . All the parameters of the { M1 } branch can be identified using manufacturer’ s data sheet. Parameters of the branches M2 and M3 though, are more difficult to obtain since they must be obtained experimentally . Furthermore, these parameters are not univ ersal and they depend on the time scale of the phenomena to be observed (described by the RC time constant τ RC = RC ). The time scales are arbitrary , howe ver they usually imply a range from sev eral minutes to several weeks or ev en more [15], [17]. B. Simplification of the supercapacitor cell model This section shows that the SC model for power system dynamics studies can be simplified and described by only the { M1 } branch and with at least one parallel group R s k C s k . The experimental data from [15], [17] are used to dev elop dedicated models in Matlab/Simulink and then use them in simulations for model simplification. The simulation results are related to two commercial SCs (Epcos 110 F and Maxwell 140) represented by different models with varying le vels of detail to validate model simplification. 3 Different model responses are produced in MA TLAB- Simulink using Simscape Electrical toolbox. The number of branches in { M2 } and { M3 } groups is being sequen- tially reduced and the different model responses to the charge/dischar ge test are compared. Input to the model is the current i sc ( t ) and output of the model is the SC terminal voltage u sc ( t ) . Results are shown in Fig. 2. For clarification, 6 branch model represents the total number of branches (first branch { M1 } , 4 parallel branches { M2 } and a self-discharge branch { M3 } ). Results for both SCs with dif ferent levels of detail (Fig. 2b and Fig. 2c) sho w that the branches { M2 } and { M3 } do not ha ve an impact on the model accuracy for the time scale of interest. Therefore, all the branches except the first branch can be neglected. In the next step, the adequate number of parallel RC groups in the first branch is determined. Results obtained for Maxwell 140 model are sho wn in Fig 3a. The significance of parallel RC groups depends on the dif ference between R dc and R s . If those resistances are very close together (e.g. < 1 m Ω as in Epcos model) then the impact on the accurac y is negligible. Ho wever , if the difference between R dc and R s is larger (e.g. > 2 m Ω as in Maxwell 140 model) then the accuracy is significantly impacted which is important from both the a vailable energy perspectiv e and control perspective. Mean absolute percentage error for the Maxwell model with different numbers of groups is shown in T able I. Since R dc and R s in reality depend on the cell in question, it can be conservati vely concluded that at least one parallel group should be included: both R s and the parallel RC groups will cause voltage drops and ener gy losses in the circuit, thus should be included. Since the only voltage that can be measured is the one across the supercapacitor terminals ( u sc ), R s and parallel RC groups will impact both the estimation of SoC and the DC voltage dynamics which determine the operating range of the supercapacitor module. Fig. 3b compares voltage response between the ideal ca- pacitor and the nonlinear model with 1 RC group. The closest voltage profile w as obtained when the capacitance of the ideal capacitor was set to the SC capacitance at half the rated voltage. Nev ertheless, the ideal representation will not reflect the voltage transient effect due to the ESR which occurs when the charging or discharging current is discontinued. This voltage transient is important from the control perspective because it impacts the logic that enables or disables the charg- ing/discharging based on state-of-voltage. More importantly , Fig. 3c shows the difference between stored energy in a nonlinear model and in the ideal capacitor . Using capacitance at rated v oltage is overly optimistic re garding the stored ener gy for the same applied voltage, while using capacitance at 0 V is closest to the nonlinear model (mean absolute percentage error in energy for the operating points between 0.5 V and 2.5 V is 27%). This ov erestimation of the charge of the ideal ultracapacitor is coming from i) Dif ferent initial stored energy in the steady-state (capacitance of the ideal capacitor is constant while real capacitor has constant and variable, i.e., proportional to voltage, component); ii) Different rate of discharge because of the variable capacitance; iii) Internal losses due to equiv alent series resistance and the parallel RC groups of the first branch in Fig. 1. (a) (b) (c) Fig. 2. Comparison of model response for dif ferent number of branches: a) T est current 30 A, b) Maxwell BCAP0140 voltage; c) Epcos 110 F voltage T ABLE I A C C U RA CY O F T H E S U PE R C A P A CI T O R M OD E L ( M A X WE L L ) W I T H D I FFE R E N T N UM B E R S O F P A R A LL E L R C G RO U P S I N T H E FI RS T B R A N CH Number of RC groups Mean absolute percentage error [%] V oltage Energy 1 (relative to 5 group model) 2.6 6.5 0 (relative to 5 group model) 6.7 12.9 0 (relative to 1 group model) 4.3 10.3 C. Building a super capacitor bank model In section II-B it was shown that the SC dynamics can be accurately represented using first branch only { M1 } with at least one parallel RC group. T o build a capacitor bank of a higher power rating, a certain number of cells n s can be connected in series to form a string and a certain number of strings n p can be connected in parallel to form a module. Modules can then be connected in parallel to form a bank. Assuming completely identical cells, it is easily shown using Kirchoff ’ s voltage and current laws that the voltage of the string u s sc and the current of the module i m sc are equal to (5) and (6), respectiv ely . u s sc ( t ) = n s u sc ( t ) (5) i m sc ( t ) = n p i sc ( t ) (6) Assumption of identical and balanced cells is reasonable for bulk power system simulations, otherwise the modelling would ha ve to be done on an electronic component level. Commercial solutions always have some sort of cell balancing system implemented to keep the cell usage and aging uniform. Cell aging is reflected in the change of parameters so that can 4 5 RC groups 1 RC group 0 RC groups (a) (b) Stored energy [J] (c) Fig. 3. Maxwel BCAP0140 model response: a) V oltage response to different numbers of RC groups; b) V oltage response compared to ideal capacitor; c) Energy stored in a supercapacitor with respect to voltage be accounted for if the parameters of an older system are known. Finally , the dynamic model of the bank can be built using circuit analysis in the time domain for the first branch only by setting u sc ( t ) as an output y ( t ) , i sc ( t ) as an input u ( t ) . Capacitor voltages are chosen as state v ariables. Complete nonlinear model of the SC bank in the analytic form is described by (7)–(12) and shown in Fig. 4 where R s k and C s k are defined by (2) and (3), respectiv ely . u sc ( t ) = i sc ( t ) R s + u C ( t ) + n X k =1 u C s k = y ( t ) (7) i sc ( t ) = u ( t ) (8) u s sc ( t ) = n s u sc ( t ) = n s y ( t ) (9) i m sc ( t ) = n p i sc ( t ) = n p u ( t ) (10) du C dt = i sc ( t ) C 0 + k v u C ( t ) (11) du C s k dt = − u C s k R s k C s k + i sc ( t ) C s k (12) Note that (11)–(12) assumes the relation dQ = C sc · du C per [17] instead of Q = C sc · u C . This is an alternativ e definition of capacitance for voltage-dependent capacitors [19]. Never - theless, in both definitions the expression for capacitance in the denominator of (11)–(12) remains linearly proportional to voltage only dif fering in the magnitude of k v . The behavior remains the same and the capacitance is varied by varying k v in the case studies presented here. SC cells are very low voltage devices (rated voltage up to 2.5–3 V) and they hav e to be scaled up to form a bank of grid-scale power rating, e.g. 1 MW to 100 MW . The number of cells necessary to form a grid-scale SC bank is in the order of 10 3 for 1 MW bank up to 10 5 for 100 MW scale. I I I . S U P E R C A P AC I T O R C O N T R O L S Y S T E M The complete control system consists of in verter PQ con- trol, charge/dischar ge control, DC current calculation and frequency control loop. The block diagram of the complete SC bank energy storage system is shown in Fig. 5. P and Q are active and reacti ve power injected or absorbed by the in verter to or from the grid, while asterisk (*) denotes a set- point v alue. V g rid ac is the A C v oltage of the bus the in verter is connected to. i d and i q are the direct and quadrature axis currents of the in verter . In verter is controlled in the grid v oltage reference frame in which the grid voltage phase is estimated by the PLL (model version 3 in PowerF actory). PLL also estimates the system frequency for frequency control block. DC current calculation block calculates the SC current for charging or dischar ging. It should be noted that the presented SC control system is similar to a battery control system [20] since the requirements are the same (constant power output). The difference from the other control schemes howe ver , is that the voltage measurement is directly used as a measure of en- ergy (State-of-V oltage, SoV) rather than State-of-Char ge (SoC) through current integration since the energy of a capacitor is directly proportional to the voltage. Furthermore, capacitors are much more sensitive to applied voltage, which varies significantly more than in batteries. Therefore, a special care must be taken not to overchar ge the SC since even the v oltage which is only 5% above the rated voltage can damage the cell. Similarly the SC should not overly discharge because of current limitation for constant power . The individual controller blocks are further elaborated in the following subsections. A. Charg e contr ol Fig. 6 shows the structure of this block. State-of-V oltage (SoV) measurement is used to control the charging and discharging process. Charging is stopped if the SC bank is charged to nominal voltage, while discharging is stopped when the SC voltage f alls belo w a user defined low voltage threshold. Charging/dischar ging is enabled again when the voltage reaches a user defined minimum voltage le vel for charging/dischar ging. The input to the block are the d and q axis currents i 0 d and i 0 q from the PQ control, while the final in verter current set-points i ∗ d and i ∗ q are determined by this block. Simple low-v oltage ride through logic and current limitation are also implemented in this block. They are not shown here since they are av ailable in many papers, e.g. [6], [20], [21]. 5 i m sc ( t ) 1 /n p R s ÷ C 0 + u C ( t ) k v i sc ( t ) 3( R dc − R s ) π 2 1 1 2 ÷ − 2 R 1 2 2 ÷ − R 1 n 2 ÷ − R P u R s ( t ) R u C ( t ) u C s 1 ( t ) u C s 2 ( t ) u C s n ( t ) n s u s sc ( t ) u sc ( t ) Fig. 4. Block diagram of the nonlinear supercapacitor module model Grid frequency control PQ control P ∗ ; P Q ∗ ; Q f i ∗ d V ∗ ; V grid ac δP Charge control & L VR T i 0 d i 0 q DC current calculation Ultracapacitor/ supercapacitor model i m sc u sc ; u s sc Inv erter i ∗ d i ∗ q i max ch i max dch P PLL cos θ ∗ ; sin θ ∗ V grid ac f P Q Fig. 5. Complete model of the supercapacitor bank system u sc ≥ U max ch disable charging flag ≤ U start ch enable charging flag ≤ U min dch disable discharging flag ≥ U start dch enable discharging flag SR flip-flops S R S R Q − ¯ Q i max ch ∈ { 0 , 1 } i max dch ∈ {− 1 , 0 } V grid ac i 0 d i 0 q i max ch i max dch L VR T logic and current limitation i ∗ d i ∗ q Fig. 6. Charge control, L VR T and current limitation block B. DC curr ent calculation The input to the SC model is the current, ho wever in po wer system applications the power is usually controlled and not the current. This block calculates the charging or discharging DC current based on the actual inv erter power output. Block diagram of this subsystem is shown in Fig. 7. It should be noted that this module as well as the SC model works with SI units, while other subsystems work in p.u. I max ch and I max dch are the maximum single cell charging and discharging current in A (e.g. ± 100 A). C. PQ contr ol Fig. 8 shows the PQ control structure of the SC bank in verter . In this case, the inv erter is modelled as a controlled current source and the d and q axis currents are obtained from the acti ve and reactive power control error , respecti vely . Mea- surement/control lag is also included in this block diagram. The term i ∗ d − i 0 d is a compensation term for acti ve po wer during low-v oltage ride through when the active power should be low and reactiv e po wer high. Reactive power or terminal voltage control can be both chosen. Howe ver , if reactiv e power control is chosen it will be ov erridden by terminal voltage control P u s sc ÷ i max ch I max ch n p i max dch I max dch n p 1 τ c s +1 Control and measurement lag i m sc Fig. 7. DC current calculation block during low-v oltage ride through. Simplified structure of the con verter and its control applicable for electromechanical transient simulations of large networks (integration step size 1–10 ms) which neglects fast inner current control loops and A C-side filter/grid dynamics was used in this study [21]. Other or more detailed models could be used for dif ferent studies, e.g. EMT , single-machine infinite bus or smaller systems [7], [9], [12]. D. Grid fr equency control This block is shown in Fig. 9. The input to this block is the grid frequency signal estimated by the PLL and the output is the requested change in po wer . The type of implemented algo- rithm for frequency response can be arbitrary . Ho wev er , based on the SC characteristics, in this study two control loops are used. The bottom loop is a standard virtual inertial response with a washout filter to make the output signal smoother since the time deriv ative operation inherently amplifies noise. The upper loop is more akin to a standard droop control, but it also has a washout filter which means this contribution will diminish in steady-state, hence the name quasi-dr oop . The reasoning for this choice is the follo wing: the SC does not ha ve a lot of stored ener gy—if the standard droop control is employed then the SC output power is initially proportional to the frequency de viation. Howe ver , once the SC is discharged, the output power will fall to zero which will cause a bigger secondary frequency drop. By setting a large washout filter time constant, the output power will slowly diminish while the con ventional units pick up. Therefore, the difference between the inertia control loop and quasi-droop control loop is in the washout filter time constant ( τ d w τ i w ). I V . S I M U L A T I O N A N D R E S U L T S The performance of the proposed model is implemented and tested on a standard IEEE 14-bus test system included in DIgSILENT PowerF actory . T wo scenarios are tested: loss of generating unit and low-v oltage ride through event. Results obtained with the ideal model and the proposed nonlinear model are compared for different capacitor sizes and initial conditions. 100 MW SC bank is connected to b us 06. Param- eters of the SC bank system are given in the Appendix. An actual supercapacitor system was not av ailable for validation hence experimentally v alidated models were used for v alidation of model reduction. There is no complete control system model of SC in the open literature that can be easily reproduced in system studies. The model proposed here is modular and generic enough to be implemented in standard power system simulation software or modified if necessary . It can adequately capture the relev ant SC dynamics as the modelled SC control is generic/flexible with the capacitor 6 P P ∗ δ P − + − 1 τ c s +1 Control and measurement lag K d p + K d i s i 0 d − i ∗ d Q Q ∗ V V ∗ − − 1 τ c s +1 K q p + K q i s i 0 q Low-voltage flag Fig. 8. Supercapacitor bank inverter PQ control Quasi-droop control Inertia con trol ∆ f K d Gain sτ d w sτ d w + 1 W ashout filter K i sτ i w sτ i w + 1 + δP Fig. 9. Supercapacitor bank grid frequency control module model based on past experimentally validated results [15], [17]. A. Loss of gener ating unit T wo cases are analysed to trigger an underfrequenc y ev ent: i) at t = 1 s, generator G02 deliv ering 95 MW is disconnected; ii) at t = 1 s, generator G02 deliv ering 190 MW is discon- nected. The performance of the ideal model is compared to the performance of the nonlinear model. There are se veral param- eters which may affect the performance of both models: initial capacitor voltage, ideal capacitor capacitance and how much the SC capacitance v aries with voltage (i.e. ratio of capacitance at 0 V to capacitance at rated v oltage). SC bank size is 100 MW . Characteristics of the system frequency response were observed and results for VIR control and quasi-droop control for 95 MW disturbance are shown separately in Fig. 10 and Fig. 11, respectively (shown separately because parameters of quasi-droop control impact the frequency response differently; nonetheless, the two control loops can be used simultaneously as shown in Fig. 9). There is no significant dif ference in frequenc y nadir between using ideal and nonlinear model for a range of initial SC voltage values and capacitance models (Fig. 10a–10c, 11a– 11b); howe ver , the actual range depends on the size of the disturbance and control system parameters. T able II shows the operating range in terms of SoC in which the ideal model does not adequately represent the nonlinear model (inaccurate range) meaning that the range in which the ideal model is indeed adequate (accurate range) is complementary to the shown range. The criterion for selection of inaccurate range is the relative error between the nonlinear and ideal model for maximum frequenc y deviation (if this error is > 5 %, then the ideal model does not adequately represent the nonlinear model). In T able II, best ideal model is highlighted by a shaded cell. Generally , for the same size of SC bank the ideal model is accurate in a wider range for a smaller disturbance, and in a more narro w range for a larger distrubance. Furthermore, smaller gain and longer washout time constant of quasi-droop control increase the accurate range of the ideal model. In all cases, the inaccurate range is wider for larger capacitance of ideal model, although the values obtained with quasi-droop control are not as sensitiv e to the value of capacitance of the ideal model. For a smaller disturbance, the inaccurate range is between 20% and 50% SoC for VIR control and between 20% and 30% for quasi-droop control meaning that the ideal model inaccurately represents the nonlinear model for low to medium SoC. When the disturbance exceeds the size of the SC, the inaccurate range is between 25% and 70% for VIR, although this range can be shifted to the higher SoC range if the v ariable part of SC capacitance is bigger (e.g. 71% to 85% SoC for 25% variable capacitance and 70% to 90% SoC for 40% variable capacitance). The inaccurate range for quasi-droop control is between 20% and 40% SoC. Note that for both disturbance sizes, there are empty cells in the table which mean that the ideal model in those cases accurately represents the nonlinear model for the whole operating range. One must k eep in mind that the results are shown for an underfrequency e vent—for an ov erfrequency e vent the observed beha viour is complementary , i.e. the inaccurate range is in higher SoC because the capacitor is char ging in that case and is limited by the maximum v oltage. Fig.10d shows a verage RoCoF for the SC with 40% variable capacitance. It can be seen that the RoCoF increases with lower SoC because the SC bank will not be able to deli ver the requested power before reaching minimum SoV . The RoCoF plots correspond to the frequency nadir plot (e.g. compare Fig. 10d and Fig. 10c) for all SC models. Generally , larger variable capacitance will result in larger RoCoF and larger nadir for the same initial conditions (e.g. compare the nadir at 1.7 V in Fig. 10). The observed differences in nadir and RoCoF between different SC expressions are < 0 . 05 Hz and < 0 . 05 Hz/s, respectiv ely . Maximum possible difference in frequency nadir for all analysed scenarios is sho wn in T able III. In 34 out of 36 analysed cases, the ideal model with minimum capacitance will yield the best results (error is ≤ 0 . 1 Hz) for both types of control. Howe ver , if the v ariable capacitance is larger (e.g. 40% model), than the smallest maximum difference was achiev ed for the ideal model with av erage capacitance ( C @0 . 5 U r ). Fig. 11c and Fig. 11d sho w the time to dischar ge of 10% and 40% variable capacitance models. The size of variable capacitance impacts which ideal model adequately describes the nonlinear model in terms of time to discharge. The mean absolute relati ve error in discharge time between the ideal model and the nonlinear model with 10% variable capacitance is 16.4% for ideal model with C @ U r , 12.8% for ideal model with C @0 . 5 U r and 8.7% for ideal model with C @0 V . The same error for 40% v ariable capacitance is 24.5%, 9.9% and 10.4%, respecti vely . As the variable part of capacitance increases, so does the error of the ideal models with minimum and maximum capacitance and the average model will be the most accurate on av erage in terms of time to discharge taking into account the whole operating range. T aking into account frequency nadir, RoCoF and time to discharge, the best ideal model for SC with up to 25% variable 7 capacitance is the ideal model with minimum capacitance ( C @0 V) and this model is the same for both types of control. Howe ver , if the v ariable capacitance is larger than that, e.g. 40%, the best ideal model may actually be the model with a verage capacitance for both types of control (for VIR control this is with respect to maximum absolute difference in frequency nadir while for quasi-droop control this is with respect to time to discharge). The observed differences are important to accurately predict the performance of a SC bank during all operating points and therefore the grid frequency dynamics during large load- generation disturbances. For example, if the frequency is to be contained after the loss of a generator, using an ideal model might lead to a conclusion that the frequency will indeed be contained, while in reality , the SC will not be able to deliv er necessary power output and underfrequency load- shedding (UFLS) will be triggered (e.g., see Fig. 10a and Fig. 11a at 1.5 V if the first stage of UFLS is set at 59.2 Hz). SC can sustain an output power profile for a limited amount of time and it is necessary to accurately estimate the stored energy in order to guarantee the requested provision of system services. Choosing an inadequate capacitance value may lead to a difference in calculated frequency of over 0.1 Hz (e.g. Fig. 10c capacitance at rated voltage) as well the RoCoF (Fig. 10d). B. Discussion on power and ener gy of a super capacitor for a loss of generating unit An ideal capacitor will extract the maximum amount of energy between two v oltage lev els without any losses, ho wev er in a real supercapacitor a part of the stored energy will be dissipated as a thermal loss that is time-variant as the discharging current changes to maintain the requested power profile. Thus, the effecti ve capacitance of an ideal model would hav e to be recalculated for each operating point which is not practical and since the real supercapacitor is nonlinear and the requested power profile depends on the external factors (i.e. size of disturbance, grid dynamics and frequency controller tuning) this is also practically impossible to do both analytically and apriori. In a general case, the ideal model will behav e identically to the realistic model in terms of output power profile if the ener gy requested during the transient is less than the minimum of stored effecti ve energy in an ideal model and in realistic model: E tran < min { E ideal , E real } (13) Considering the same loss of a 95 MW generating unit G02 at t = 1 s, the power profile and SoC are observed for a fully (100% SoC) and partially ( ∼ 45% SoC) charged super- capacitor for different ideal model representations. Frequenc y controller parameters are set to K i = 100 p.u., K d = 50 p.u., τ i w = 1 s, τ d w = 30 s to induce a complete dischar ge of the supercapacitor bank. Fig. 12 shows the results for a supercapacitor with 10% variable capacitance. It can be seen that in the case of a fully charged supercapacitor , the model with av erage capacitance is the most accurate in terms of power profile and SoC (Fig. Frequency nadir [Hz] (a) C sc = 900 + 40 u C ; C 0 /C max = 0 . 9 Frequency nadir [Hz] (b) C sc = 750 + 90 u C ; C 0 /C max = 0 . 75 Frequency nadir [Hz] (c) C sc = 600 + 150 u C ; C 0 /C max = 0 . 6 Average RoCoF [Hz/s] (d) A verage RoCoF for C sc = 600 + 150 u C Fig. 10. Frequency nadir with different supercapacitor models for virtual inertial response control Frequency nadir [Hz] (a) C sc = 900 + 40 u C ; C 0 /C max = 0 . 9 Frequency nadir [Hz] (b) C sc = 600 + 150 u C ; C 0 /C max = 0 . 6 Time to discharge [s] (c) C sc = 900 + 40 u C ; C 0 /C max = 0 . 9 Time to discharge [s] (d) C sc = 600 + 150 u C ; C 0 /C max = 0 . 6 Fig. 11. Frequency nadir and time to discharge with different supercapacitor models for quasi-droop control 12a and Fig. 12b): ideal model with av erage capacitance can sustain the output power only for a fraction of a second longer than nonlinear model (14.85 s compared to 14.38, or 3.3%). On the other hand, ideal model with rated capacitance sustains the output po wer 1.6 s longer (15.99s compared to 14.38 s or 11.2% longer), while the ideal model with minimum capacitance sustains the output po wer 0.6 s shorter (13.77 s compared to 14.38 s, or -4.2%). Fig. 12b sho ws that energy depletes faster in the case of ideal model with minimum ca- pacitance and slower in the case of ideal model with maximum 8 T ABLE II O P ER ATI N G R A N GE ( S O C ) I N W H IC H I D E A L M OD E L D O E S N OT D E S CR I B E N O NL I N E AR M O DE L AC C U RAT ELY W I TH R E SP E C T T O M A XI M U M F R EQ U E NC Y D E VI ATI O N ( > ± 5 % E R RO R ) ; NA R RO WE R R A N G E I S B E TT E R Designed operating range: 1.1 V – 2.7 V (40%–100% U r ∼ 15%–100% SoC) Disturbance size: 95 MW Model VIR control Quasi-droop control C @0 V C @0 . 5 U r C @ U r C @0 V C @0 . 5 U r C @ U r C 0 /C max = 0 . 9 22%–30% 22%–43% 22%–43% 22%–30% 22%–30% 22%–30% C 0 /C max = 0 . 75 21%–28% 21%–33% 21%–42% 21%–28% 21%–28% 21%–28% C 0 /C max = 0 . 6 40%–50% 20%–35% 20%–45% - 20%–27% 16%–27% Disturbance size: 190 MW Model VIR control Quasi-droop control C @0 V C @0 . 5 U r C @ U r C @0 V C @0 . 5 U r C @ U r C 0 /C max = 0 . 9 26%–43% 26%–48% 26%–66% 26%–39% 21%–39% 21%–39% C 0 /C max = 0 . 75 71%–85% 25%–47% 25%–71% - 25%–42% 21%–42% C 0 /C max = 0 . 6 70%–90% 23%–45% 23%–70% - 23%–40% 20%–40% T ABLE III M A XI M U M A B SO L U T E D IFF E R E NC E I N F R E Q UE N C Y NA D I R F O R A LL A NA LYS E D C A S ES ( S MA L L E R I S B ET T E R ) Disturbance size: 95 MW Model VIR control Quasi-droop control C @0 V C @0 . 5 U r C @ U r C @0 V C @0 . 5 U r C @ U r C 0 /C max = 0 . 9 0.060 Hz 0.067 Hz 0.072 Hz 0.066 Hz 0.069 Hz 0.071 Hz C 0 /C max = 0 . 75 0.046 Hz 0.062 Hz 0.082 Hz 0.050 Hz 0.065 Hz 0.074 Hz C 0 /C max = 0 . 6 0.069 Hz 0.057 Hz 0.103 Hz 0.035 Hz 0.059 Hz 0.077 Hz Disturbance size: 190 MW Model VIR control Quasi-droop control C @0 V C @0 . 5 U r C @ U r C @0 V C @0 . 5 U r C @ U r C 0 /C max = 0 . 9 0.096 Hz 0.104 Hz 0.117 Hz 0.101 Hz 0.113 Hz 0.124 Hz C 0 /C max = 0 . 75 0.089 Hz 0.101 Hz 0.134 Hz 0.080 Hz 0.102 Hz 0.130 Hz C 0 /C max = 0 . 6 0.149 Hz 0.094 Hz 0.167 Hz 0.065 Hz 0.099 Hz 0.137 Hz capacitance, while the ideal model with a verage capacitance is fairly accurate in this case. When the supercapacitor is only partially charged when the disturbance happens, the differences are more pronounced (Fig. 12c and Fig. 12d) and all ideal models o verestimate the available energy: ideal models with minimum, average and rated capacitance sustain the output po wer for 0.5 s (+13.9%), 0.74 s (+20.6%) and 0.98 s (+27.3%) longer , respecti vely . Furthermore, initial SoC is overestimated by ideal models by 1.1%, it depletes slower and total change of SoC is 15.4% bigger using ideal models. Fig. 13 shows the results for a supercapacitor with 40% variable capacitance. In the case of a fully charged the most accurate ideal model in terms of power profile and SoC is again the ideal model with av erage capacitance (Fig. 13a and Fig. 13b): it can sustain the output power profile for 1 s shorter (11.88 s compared to 12.95 s or -8.3%), while the ideal model with minimum capacitance can sustain the output power for 4.5 s shorter (8.43 s compared to 12.95 s, or - 35.1%) and the ideal model with maximum capacitance can sustain the output po wer for 3 s longer (15.9 s compared to 12.95 s, or +22.8%). Fig. 13b sho ws that the SoC can deplete at significantly different rates depending on the choice of ideal capacitance. For a partially charged supercapacitor (Fig. 13c and Fig. 13d), the most accurate ideal model with respect to power profile and SoC is the one with minimum capacitance: it discharges 0.22 s earlier (-7.0%) than the nonlinear model. Ideal model with average and rated capacitance overestimate the stored energy and they can sustain the po wer output for a longer time: 0.54 s longer (+17.1%) and 1.37 s longer (+43.3%). Ideal model ov erestimates the initial SoC by 4.6% and the total change of SoC is 21.4% bigger using ideal model. It can be concluded that as the variable part of capacitance increases, there is a bigger error in estimation of initial SoC as well as in depletion rate of energy between ideal models and nonlinear model. Also, as the initial SoC decreases, different ideal models describe the nonlinear model better . Therefore, the nonlinear v oltage dynamics and losses of a real supercapacitor cell make the modelling using ideal or simplified models problematic. C. Low voltage ride through Performance of different models has also been ev aluated for balanced three-phase faults near Bus 06. In the first scenario, a 500 ms three-phase self-clearing fault has been applied to Bus 02 and Fig. 14 shows the SC DC v oltage for different initial voltages and SC models. In all cases, the initial voltage and SC voltage dependence do not influence the model performance significantly . The difference between the ideal and nonlinear model (IM,NM) is more visible at the first glance: the nonlinear model has a greater voltage change due to the voltage drop on the ESR. Ho wev er , maximum difference is less than 10 V so it can also be neglected. Bus 06 A C voltage profile is sho wn in Fig. 15. Similarly , the initial voltage, SC voltage dependence and the type of model do not have a significant impact on the post-fault voltage transient, while the profiles are identical during the fault. Maximum observ ed dif ference of absolute v oltage v alues is 0.03 p.u. between nonlinear and ideal model (both fully charged) for both SC capacitance expressions. For a partially charged SC ( U 0 = 2 V), this difference is 0 for both capacitance expressions. This is because the system is initially accelerating after the fault clears, so the SC is char ging and the fully charged SC is the worst case scenario since it can’t accept much more char ge (a small tolerance of < 1% e xists between the rated voltage and cut-of f voltage). Therefore, SC with low SoC will behav e identically to the partially charged SC (e.g., U 0 = 2 V or about 50% SoC). Maximum observed difference of absolute voltage v alues between partially charged and fully charged SC is 0.07 p.u. The impact of dif ferent faults and durations ha ve also been analyzed. Fig. 16a shows the difference in model performance 9 Power [MW] (a) Power profile for full charge SoC [%] (b) SoC profile for full charge Power [MW] (c) Power profile for partial charge SoC [%] (d) SoC profile for partial char ge Fig. 12. Power and SoC profiles for fully charged (top) and partially charged (bottom) supercapacitor for 10% variable capacitance ( C sc = 900 + 40 u C ) Power [MW] (a) Power profile for full charge SoC [%] (b) SoC profile for full charge Power [MW] (c) Power profile for partial charge SoC [%] (d) SoC profile for partial char ge Fig. 13. Power and SoC profiles for fully charged (top) and partially charged (bottom) supercapacitor for 40% variable capacitance ( C sc = 600 + 150 u C ) for dif ferent duration (100, 200 and 300 ms) of a fault resulting in a 37% voltage dip at Bus 06. The initial SoC of the SC is 76–78%. It can be seen that the nonlinear model will hav e a somewhat higher voltage spike and a more oscillatory behaviour after the fault clears due to the ESR. Practically , this means that the overv oltage protection may be triggered sooner when a nonlinear model is used compared to an ideal model. Again, the difference between the IM and NM is in tens of V and will not play a significant role in grid dynamics. The longer the fault duration, the bigger the v oltage spik e and the bigger the difference between models. The size of the voltage dip (20% to 80%) has a similar impact (Fig. 16b). Ideal model will hav e smaller and smoother transients than the nonlinear mode due to the lack of ESR, but the observed dif ferences will not significantly impact the grid dynamics. Same as before, depending on the operating point of the SC, the undervoltage or ov ervoltage protection may be triggered sooner when the nonlinear model is used. In summary , the ideal model representation of SC is ad- equate for transient stability simulations. The difference in internal dynamics behavior is small enough that it should not hav e any impact on the grid results. V . C O N C L U S I O N In this paper, an accurate supercapacitor bank model and associated control system has been presented for the use in power system dynamics simulations. Starting from the most detailed RC model of a supercapacitor cell, the model has been gradually reduced until arriving to the simplest representation which adequately describes the supercapacitor dynamics, confirmed by simulation experiments. The proposed model is described with only 4 parameters which are easily obtained from manufacturer’ s data sheet: capacitance at zero voltage, voltage-dependent capacitance, DC resistance and high-frequency resistance. The performance of the presented model compared to an ideal model has been tested in an IEEE 14-bus test system in frequency control and L VR T scenarios. For frequency control, the ideal model does not always represent the nonlinear model adequately depending on the initial supercapacitor v oltage and disturbance size. For an un- derfrequency ev ent, a fully to partially charged supercapacitor may be adequately represented by an ideal model in terms of system frequency response, but nearing the minimum voltage limit the ideal model may yield o verly optimistic or pessimistic results (frequency nadir difference of ov er 0.1 Hz can be observed depending on the ideal capacitor capacitance v alue). Similar behavior is observed for an ov erfrequency event. Generally , equiv alent series resistance and parallel RC groups of the first branch reduce the ef ficiency of the supercapacitor , while voltage-dependent capacitance changes the amount of stored energy during charging/dischar ging and influences the charge/dischar ge rate. The observed mean absolute relati ve error in discharge time between ideal and nonlinear model ranges from 9% to 16% for a 10% variable capacitance and between 10% and 25% for 40% v ariable capacitance, while the maximum observed absolute relativ e error in discharge time can go up to 27% for 10% v ariable capacitance and 43% for 40% variable capacitance. The best ideal model for most observed cases for both types of presented control schemes is the ideal model with capacitance set between minimum and average supercapacitor capacitance, with the ideal model with average capacitance being usually more accurate for high initial SoC. Losses and nonlinear voltage dynamics complicate representing a real supercapacitor with an ideal model for all operating points. 10 2.6 2.8 3 3.2 Time [s] 990 995 1000 720 725 730 735 740 DC voltage [V] (a) C sc = 900 + 40 u C ; C ideal = 900 F DC voltage [V] (b) C sc = 600 + 150 u C ; C ideal = 600 F Fig. 14. DC voltage of a supercapacitor bank during fault for different initial voltages and supercapacitor models Bus 06 voltage [p.u.] (a) C sc = 900 + 40 u C ; C ideal = 900 F Bus 06 voltage [p.u.] (b) C sc = 600 + 150 u C ; C ideal = 600 F Fig. 15. A C voltage at Bus 06 after a disturbance DC voltage [V] (a) C sc = 600 + 150 u C ; C ideal = 600 F DC voltage [V] (b) C sc = 600 + 150 u C ; C ideal = 600 F Fig. 16. DC voltage of a SC bank for a): different fault duration τ f for a 37% voltage dip; b): different sizes of the voltage dip ∆ u f for a 200 ms fault For low-v oltage ride through, the impact of modelling is not significant and the ideal model will be adequate, although the undervoltage and overv oltage protectiv e circuits may be triggered sooner for the nonlinear model. A P P E N D I X n s = 370, n p = 400 , bank rated power: 100 MW , C max ≈ 1000 F , R dc = 0 . 5 m Ω , R s = 0 . 25 m Ω , I max ch /I max dch = ± 615 A , U max ch = 2 . 71 V , U start ch = 2 . 4 V , U min dch = 1 . 1 V , U start dch = 1 . 4 V , τ c = 50 ms, K i = 150 p.u., K d = 100 p.u., τ i w = 1 s, τ d w = 30 s, K d p = K q p = 1 p.u., K d i = K q i = 100 p.u. R E F E R E N C E S [1] X. Luo, J. W ang et al. , “Overview of current development in electrical energy storage technologies and the application potential in power system operation, ” Applied Energy , vol. 137, pp. 511–536, jan 2015. [2] Maxwell T echnologies, “Ultracapacitors frequency response application brief, ” 2018, technical brochure. [3] A. T ahri, H. E. Fadil et al. , “Management of fuel cell power and supercapacitor state-of-charge for electric vehicles, ” Electr . P ower Syst. Res. , vol. 160, pp. 89 – 98, 2018. [4] T . Zhou and W . Sun, “Optimization of battery–supercapacitor hybrid energy storage station in wind/solar generation system, ” IEEE T rans. Sustain. Energy , vol. 5, no. 2, pp. 408–415, April 2014. [5] M. F . M. Arani and E. F . El-Saadan y, “Implementing virtual inertia in dfig-based wind po wer generation, ” IEEE T rans. P ower Syst. , vol. 28, no. 2, pp. 1373–1384, May 2013. [6] S. I. Gkavanoudis and C. S. Demoulias, “ A combined fault ride-through and power smoothing control method for full-con verter wind turbines employing supercapacitor energy storage system, ” Electr . P ower Syst. Res. , vol. 106, pp. 62 – 72, 2014. [7] J. F ang, Y . T ang et al. , “ A battery/ultracapacitor hybrid energy storage system for implementing the power management of virtual synchronous generators, ” IEEE T rans. P ower Electron. , v ol. 33, no. 4, pp. 2820–2824, April 2018. [8] M. H. Fini and M. E. H. Golshan, “Determining optimal virtual inertia and frequency control parameters to preserv e the frequency stability in islanded microgrids with high penetration of renew ables, ” Electr . P ower Syst. Res. , vol. 154, pp. 13 – 22, 2018. [9] V . Gevor gian, E. Muljadi et al. , “Supercapacitor to provide ancillary services, ” in 2017 IEEE Energy Con version Congress and Exposition (ECCE) , Oct 2017, pp. 1030–1036. [10] L. Sigrist, I. Egido et al. , “Sizing and controller setting of ultracapacitors for frequenc y stability enhancement of small isolated power systems, ” IEEE Tr ans. P ower Syst. , vol. 30, no. 4, pp. 2130–2138, July 2015. [11] J. Cao, W . Du et al. , “Optimal sizing and control strate gies for hybrid storage system as limited by grid frequency deviations, ” IEEE Tr ans. P ower Syst. , vol. 33, no. 5, pp. 5486–5495, Sep. 2018. [12] M. G. Molina and P . E. Mercado, “Modeling of a DST A TCOM with ultra-capacitor energy storage for po wer distribution system applica- tions, ” in XIII Eriac D ´ ecimo T ercer Encuentro Regional Iber oamericano de CIGR ´ E , May 2009, pp. 1–8. [13] G. Delille, B. Francois, and G. Malarange, “Dynamic frequency control support by energy storage to reduce the impact of wind and solar generation on isolated power system’ s inertia, ” IEEE T ransactions on Sustainable Energy , vol. 3, no. 4, pp. 931–939, 2012. [14] I. Egido, L. Sigrist et al. , “ An ultra-capacitor for frequency stability enhancement in small-isolated po wer systems: Models, simulation and field tests, ” Applied Energy , vol. 137, pp. 670 – 676, 2015. [15] V . Musolino, L. Piegari, and E. T ironi, “New full-frequency-range supercapacitor model with easy identification procedure, ” IEEE Tr ans. Ind. Electron. , vol. 60, no. 1, pp. 112–120, Jan 2013. [16] F . Rafik, H. Gualous et al. , “Frequency , thermal and voltage superca- pacitor characterization and modeling, ” Journal of P ower Sources , vol. 165, no. 2, pp. 928 – 934, 2007, iBA – HBC 2006. [17] R. Faranda, M. Gallina, and D. T . Son, “ A new simplified model of double-layer capacitors, ” in 2007 International Conference on Clean Electrical P ower , May 2007, pp. 706–710. [18] S. Buller, E. Karden et al. , “Modeling the dynamic behavior of superca- pacitors using impedance spectroscopy , ” IEEE T rans. Ind. Appl. , vol. 38, no. 6, pp. 1622–1626, Nov 2002. [19] I. Zeltser and S. Ben-Y aakov, “On spice simulation of voltage-dependent capacitors, ” IEEE T rans. P ower Electr on. , vol. 33, no. 5, pp. 3703–3710, 2018. [20] DIgSILENT GmbH, “Battery energy storing system template, ” 2017, template documentation. [21] WECC REMTF, “WECC wind plant dynamic modeling guidelines, ” WECC, T ech. Rep., April 2014. Matej Krpan (S’17) received the B.Eng. and M.Eng. degrees in electrical engineering from University of Zagreb in 2014 and 2016, respectively . He is currently pursuing the PhD degree at the Uni versity of Zagreb F aculty of electrical engineering and computing, Zagreb, Croatia. Igor Kuzle (S’94–M’97–SM’04) is a Full Professor at the Department of Energy and Power Systems of the University of Zagreb Faculty of electrical engineering and computing, Zagreb, Croatia. Ana Radovanovi ´ c (S’18) received the B.Sc. and M.Sc. degrees in electrical engineering and computer science from the Univ ersity of Belgrade, Belgrade, Serbia, in 2013 and 2014, respectively . She is currently pursuing the Ph.D. degree at the Department of Electrical and Electronic Engineering, the Univ ersity of Manchester, Manchester, U.K. Jovica V . Milanovi ´ c (M’95–SM’98–F’10) is currently a Professor of Electri- cal Po wer Engineering and a Deputy Head of Department of Electrical and Electronic Engineering at the University of Manchester , Manchester , U.K.
Original Paper
Loading high-quality paper...
Comments & Academic Discussion
Loading comments...
Leave a Comment