Transient behaviour veriﬁcation and controller tuning for an uncertain grid-connected photovoltaic system using reachability analysis

Power electronics–based converters for photovoltaic (PV) systems are susceptible to over-currents; it is important to design their controllers to reduce the transient current for all viable operating conditions. To design a current controller and ﬁnd the maximum transient current via simulation-based techniques, the exact values of the system parameters, initial states, and inputs are required. However, they are not precisely known in practice, some system parameters such as inductances may change over time, and output power and load are variable. The uncertainty in the parameter (ﬁlter inductance) and input of the system (injected power) should be considered in the analysis of a PV system controllers as it can degrade their performance, which are designed for the system nominal parameters. This paper employs reachability analysis for a grid-connected PV system to (1) ﬁnd the maximum transient current, (2) devise an improved PI current controller and (3) compare the maximum transient current in PI- and internal model control (IMC)-based controllers with uncertain-but-bounded input power and inductance error. Simulation and experimental studies showcase the results.


INTRODUCTION
Renewable energy resources (solar and wind) generate about 8.63% (in 2019) of the electricity in the United States [1]. Solar photovoltaic (PV) energy is a popular renewable resource, and its steady-state and transient performance analysis are of great interest. A PV unit is commonly connected to the power grid via a voltage-sourced converter (VSC) and an RL filter. To control the current of a grid-connected PV unit, a local current controller should be used. The role of the local controller is to effectively follow its set point by correctly applying the actuation variables (e.g. voltage) to the VSC. To study the behaviour of a PV system and design its controller, the exact model of the plant, its operating condition, and its inputs are needed; however, the plant parameters are not necessarily constant and may change during the system operation which may degrade the controller performance. Moreover, the injected power of the PV system depends on the sunlight which is variable; the load is also uncertain. Collectively, this uncertainty in system parame-ters and inputs introduces two important challenges related to the dynamic behaviour of a PV system: (1) design of the controller architecture and parameters; and (2) analysis of its transient behaviour. This paper reviews and addresses both of these challenges. First, to design the controller for an uncertain system, robust control methods are conventionally used. Ref. [2] proposes a robust online PID controller tuning method to update its parameters while avoiding instability during a change in the system parameters. Since this method continuously retunes the controller parameters based on the variation in system parameters, any error in measurements may cause an inaccurate tuning. Ref. [3][4][5] design robust controllers for VSC-based systems but they have complicated control procedures with many parameters to be tuned. Ref. [3] designs a controller to improve power quality, decrease harmonics, and increase stability margin during disturbances and uncertainties; however, it does not provide a proper transient performance considering a parameter uncertainty and it has nine tuned parameters. Therefore, an improved and simple controller design can be a preferred method.
Second, to accurately analyse the behaviour of a PV system (this paper focuses on the transient current behaviour), all of its parameters, initial states, and inputs are needed. The conventional method to analyse a PV system transient behaviour is to perform several simulation case studies [6][7][8]. However, gaining confidence into the performance of an uncertain system is difficult, and it is not possible to simulate infinitely many possible uncertain scenarios.
To address these shortcomings, this paper uses reachability analysis (RA) to analyse and tune the controller parameters of an uncertain PV system and improve its transient behaviour. Reachability analysis is a computational method to calculate the reachable set (operating boundaries) of the system variables including all possible scenarios for uncertain-but-bounded parameters [9,10]. This analysis can be used to ensure that the system variables do not exceed their limits. Several research efforts have employed RA to analyse uncertain power electronic systems. Refs. [10] and [11] use RA to verify voltage ride-through capability of PMSG and DFIG wind turbines, respectively. In these references, the reachable sets of the system variables are found considering grid voltage disturbances. Ref. [12] studies the droop control performance in DC microgrids using RA. In ref. [13], the large signal behaviour of a DC-DC boost converter is verified by RA. Ref. [14] uses RA to compare two different closed-loop controllers for a buck converter. In ref. [15], the transient behaviour of a DC-DC converter using hybrid automaton modelling is verified by RA, and its results are compared with Monte Carlo simulations. Monte Carlo-based simulations are used in systems with uncertainty using a large set of random scenarios [15,16]. Ref. [15] concludes that RA has faster and more confident results than Monte Carlo-based simulations. To sum up, refs. [10][11][12][13][14] use RA for transient analysis of power electronic systems with uncertain inputs, and [15] uses RA for analysis of systems with uncertain parameters for reliability purposes.
To the best of the authors' knowledge, no research paper has employed reachability analysis for a PV system for the purpose of tuning its controller parameters and analysis of its transient performance with both uncertain inputs and parameters. In the present paper, reachability analysis is employed to find the maximum transient value of the VSC terminal current considering uncertain-but-bounded DC power (input) and filter inductance (parameter) error (due to, e.g. saturation, aging, and heat). In this study, three different current controllers are considered: (1) standard PI current controller designed based on conventional methods such as Bode plot method for nominal system parameters, which does not consider system uncertainty; (2) improved PI current controller designed for the system with inductance error using our proposed RA-based method, which considers system uncertainty in controller design procedure; and (3) internal model control (IMC)-based controller [17,18]. The contributions of this paper are: • Derive the state-space model of the PV system with PI and IMC controllers considering uncertain filter inductance; • Design a method to retune the PI current controller to lower the maximum transient current in presence of inductance error; • Analyse the PV system with accurate and uncertain inductance for three different controllers using RA.
The RA results verify that in the presence of inductance error, improved PI controller designed by RA and IMC controller lead to a lower maximum transient current than the conventional PI controller. Moreover, simulation results show that they have lower transient current compared with a robust method proposed in [3]. Therefore, the system with improved PI and IMC controllers can better meet the system requirements.
This paper is organised as follows. Section 2 introduces RA. Section 3 discusses the uncertain grid-connected PV unit statespace model and controller design procedure. Section 4 presents the RA results. Section 5 validates the RA results by simulation and experimental results. Section 6 concludes the paper.

REACHABILITY ANALYSIS OF A LINEAR SYSTEM
To perform reachability analysis for a linear dynamical system, we start with its state-space model: where A is the uncertain state matrix of the system (A ∈  ⊂ R n×n ), x(t ) is the state vector, and v(t ) is the input vector (v(t ) ∈  ⊂ R n ). The initial state vector of the system is x 0 ∈  0 ⊂ R n . The input and initial state vectors are uncertain but bounded. The solution of (1) for all possible trajectories gives the reachable set (operating boundaries) of the system variables. The exact reachable set  e (t = t f ) of a dynamic system is the solution of (1) from time 0 to t f (t f is the final time when the system reaches the steady state); however, calculation of the exact reachable set is not straightforward and the overapproximated reachable set (t ) is calculated instead ( e (t ) ⊆ (t )) [10]. The union of all overapproximated reachable sets at time t k (t k = k , > 0, k = 1, … , Several methods can be used to compute the overapproximated reachable sets. Examples are polytopes [19], griddy polyhedra [20], ellipsoids [21], oriented rectangular hulls [22], zonotopes [10,11,23], support functions [24], and level sets [25]. In this paper, the zonotope-based method is used due to its closedness under linear transformation and Minkowski sum [23]. A zonotope is a polytope formed by the Minkowski sum of a zonohedron's line segments in any dimension and is defined as where g (i ) s are generators of the zonotope and c is the center of the zonotope (c and g (i ) s are vectors in R n ). Zonotope  is denoted as  = (c, ⟨g (1) , g (2) , … , g (p) ⟩).
The Minkowski sum of zonotopes  1 and  2 is defined as 1 , … , g 2 , g 2 , … , g (5) Using the superposition principle, the reachable set of the linear continuous-time system (1) can be calculated in two steps. First, the reachable set of the system with zero input and nonzero initial state is calculated (homogeneous solution  h (t )). Then the reachable set of the system considering non-zero input and zero initial state is calculated (inhomogeneous solution  f (t )). The homogeneous and inhomogeneous solutions are summed in each time step to obtain the reachable set for each time interval: The final reachable set in [0, t f ] is the union of reachable sets in each time interval and can be calculated from (2). Figure 1 shows the algorithm for reachability analysis.

Homogeneous solution
The homogeneous solution of a linear time-invariant system in the first time step is  h (t = ) = e   0 , where  0 is the set of possible initial states [9]. The term e At is calculated using Taylor series as where is the number of Taylor terms. E ( ) is added due to overapproximation error and is defined as where = and is chosen such that < 1. [−1, 1] is an interval matrix with its elements in [−1, 1] [9]. The reachable set Set up initial states x 0 , state matrix A, and input matrix v(t) from state-space model of the system Start Obtain zonotope-based model of the system based on its state-space Calculate the initial reachable set from (10) Calculate the homogeneous and inhomogeneous reachable set from (11) and (14) Calculate the reachable set in each time step by adding homogeneous and inhomogeneous sets from (6) Calculate the final reachable set as the union of the reachable sets of each time step. Set up initial states x 0 , state matrix A, and input matrix v(t) from state-space model of the system Start Obtain zonotope-based model of the system based on its state-space Calculate the initial reachable set from (10) Calculate the homogeneous and inhomogeneous reachable set from (11) and (14) Calculate the reachable set in each time step by adding homogeneous and inhomogeneous sets from (6) Calculate the final reachable set as the union of the reachable sets of each time step.   Figure 2.
The convex hull is enlarged to contain all reachable sets. This enlargement is defined as a correction matrix  [9]: The homogeneous reachable sets of the system for the first and subsequent time steps are where CH denotes the convex hull. Figure 2 shows the computation of the homogeneous reachable set in the first time step.

Inhomogeneous solution
To find the inhomogeneous solution of the system, it is assumed that v(t ) ∈  and r = − t . Therefore, The input of the system is uncertain. Therefore, an overapproximation should be performed to find the inhomogeneous reachable set of the system [9]: where E ( ) is added for overapproximation [9]. At each time step, the reachable set due to the inputs is calculated as Figure 3 shows the schematic diagram of a grid-connected PV system. The PV unit and the load are connected to the grid via a VSC and an RL filter (which may also represent a transformer). The conventional control method for VSC uses a nested-loop structure; the inner loop uses a fast PI controller to control the current, and the outer loop uses a slower PI controller to control the voltage of the DC capacitor V c by setting the d -axis reference current [26]. In this paper, the q-axis reference current is chosen randomly within the limits calculated based on the d -axis reference current. As an alternative for the inner loop, the IMC method can be used to control the currents while V c

STUDY SYSTEM DYNAMIC MODELLING AND CONTROLLER DESIGN
Current controller of the IMC-based VSC system is controlled by a PI controller [17]. Figures 4 and 5 show the block diagram of the PI and IMC controllers, respectively. In this paper, the performance of both current controllers is compared by reachability analysis. Moreover, a PI controller design methodology based on RA is proposed to reduce the maximum transient current when the inductance value is uncertain. The filter inductance value can have an error due to, e.g. aging, heat, and saturation, which should be considered in the system analysis and controller design. The filter resistance is the equivalent series resistance of the inductor which is relatively constant. To apply RA, the state space model of a system is needed, which is obtained next. Figure 4 shows the block diagram of the current controller of the VSC with a PI controller. The dynamic equations of the VSC in the dq-frame are [26,27] L f d dt

State-Space model of the VSC with PI controller
where i t,d and i t,q are the VSC terminal currents in the dq-frame, v t,d and v t,q are the terminal voltages, v s,d and v s,q are the grid voltages, u d and u q are controller commands, R f is resistance of the filter, and L f is the actual inductance of the filter. Equation (15) shows that the dq-frame currents dynamic equations are coupled due to L f terms. To decouple them, u d and u q are defined as two new control inputs where L f ,nom is the nominal inductance of the filter (possibly different from actual value of inductance L f ). Therefore, the new set of equations for dynamic model of the VSC system with inductance error can be derived as ] .
(17) Errors e d and e q are fed to a PI controller. Two internal state variables z d and z q are defined as integrator outputs. The controller commands can be calculated as where k p is the proportional coefficient of the PI controllers. By substituting e d = i d,re f − i t,d and e q = i q,re f − i t,q from Equations (18) and (19) in Equation (17), the complete state-space model of the system with PI controllers considering inductance error is written as d dt where k i is the integrator coefficient of the PI controllers. Equation (20) shows that if L f ,nom ≠ L f , small coupling terms between d -and qaxes can degrade the controller performance. i d,re f and i q,re f are the inputs of the PV system, and they are proportional to real power and reactive power, respectively (assuming v s,q = 0).

State-Space model of the VSC with IMC controller
Ref. [17] proposes a model-based alternative (IMC) to the conventional PI-based current controller for a VSC. The IMC approach considers the system model in the controller structure and has a simple structure as shown in Figure 5. Ref. [17] shows that compared with the PI-based current controllers, IMC-based current controllers have better performance with a smaller overshoot, faster step response, higher axis decoupling, and higher robustness against faults. However, [17] validates these with limited case studies. In this paper, the maximum transient current of the IMC controller is compared with that of the PI controller in all possible conditions using reachability analysis.
The state-space model of the VSC with IMC controller is needed for reachability analysis. Based on Figure 5, the IMC control inputs are defined as u ′ d , u ′ q , u ′′ d and u ′′ q : where i are the parameters of two PI controllers in IMC defined in [17]. They depend on nominal system parameters. However, the system parameters are not constant and the resulting controller parameters may not be optimal under all operating conditions. We have By substituting Equation (21) and (22) in Equation (23), the state-space model of the system with IMC controller is obtained as d dt IMC controller parameters are calculated based on a single tuning parameter [17].

PI controller design procedure using reachability analysis
Conventional VSC current controllers are PI-based, which can be designed based on well-known methods such as Bode plot, Nichols, Nyquist, and root locus. This paper uses Bode plot method for the initial controller design. These methods use nominal values of the system parameters. However, any changes in system parameters or error in their values may degrade the performance of the controller. To address this issue, this paper uses RA to retune the PI controller parameters considering an error in the filter inductance. RA, as discussed in Section 2, is performed on Equation (20)  This algorithm obtains the improved parameters in two tuning stages: (1) a coarse tuning stage to find the vicinity of the improved parameters, and (2) a fine tuning stage. In the coarse tuning stage, k p and k i parameters of the PI controller are increased from zero to k p,max and k i,max in steps of Δk p and Δk i (k p,max , k i,max , Δk p , and Δk i are heuristically chosen as 10k p0 , 10k i0 , 0.2k p0 and 0.4k i0 , respectively, where k p0 and k i0 are the initial PI controller parameters designed by Bode plot method). The reachable set is recorded for each value of k p and k i . The k p and k i values that give the lowest current limit are named k p,coarse and k i,coarse . Figure 6 shows the flowchart of the proposed PI parameter coarse tuning procedure. In the fine tuning stage, the vicinity of k p,coarse and k i,coarse , i.e. [k p,coarse − Δk p , k p,coarse + Δk p ] and [k i,coarse − Δk i , k i,coarse + Δk i ], is searched with a smaller step of 0.1Δk p and 0.1Δk i . The k p and k i values that give the smallest transient current bound are considered as the improved PI controller parameters, k p,imp and k i,imp .

REACHABILITY ANALYSIS RESULTS
In this section, reachability analysis is performed on the gridconnected PV system of Figure 3 with parameters of Table 1.
Case studies are presented for two different conditions: system with inductance error (with an error up to ±50%) and system with accurate inductance. Although error can vary based on the system conditions, the procedure of the analysis and controller design remains the same. For both case studies, RA finds the largest transient current in all possible scenarios and compares the performance of current controllers using the models in Equations (20) and (24). The system inputs i d,re f and i q,re f can take any value subject to i 2 d,re f + i 2 q,re f ≤ I 2 nom . To repre-FIGURE 6 PI controller parameter coarse tuning procedure for the system with inductance error using RA. The fine tuning procedure is similar except that Δk p and Δk i values are smaller (as mentioned in text, one-tenth of their values for coarse tuning). All variables are defined in Section 3.3. sent inputs as a zonotope, this circle is overapproximated by an octagon, since an octagon (as a zonotope) is the Minkowski sum of a zonohedron's line segments. The case studies are performed in MATLAB/CORA [28] from zero to 0.2 s with a time step of 200 and 50 s for the PI-and IMC-based controllers, respectively (IMC controller needs a smaller time step to converge).

Case II (Improved PI controller)
In this case study, to reduce the maximum transient current, the PI controller parameters are improved by RA based on the procedure discussed in Section 3.3. Figure 7(a) shows that the boundary of the transient current is 6 A, which means that reachability-based design reduces the maximum transient current boundary by 22.5% from 6.9 A.

Case III (IMC controller)
This case studies the system with IMC controller. Figure 7(a) shows that the boundary of the transient current is 5.4 A, a decrease of 37.5% from the initial PI controller and 15% from the improved PI controller. RA results show that the maximum transient current in the initial PI controller decreases with the improved PI controller, but it is not better than the IMC controller. Therefore, the improved PI controller designed by RA and IMC controller are more capable of meeting the system requirements than the initial PI controller designed based on conventional methods. Accordingly, to lower the transient current in the system with inductance error, it is a design choice to only retune the PI controllers (using improved PI controller) or to consider a different control strategy (IMC controller).

System with accurate inductance
This section discusses RA results in the system with accurate inductance. The case studies are similar to Section 4.1. By retuning PI controllers, the transient current boundary increases in the system without inductance error (the limit of the current is higher in Case B-II than in Case B-I). Therefore, with no inductance error, Bode plot method is a sufficient tool to design the PI parameters (comparable to IMC controller). However, in the presence of inductance uncertainty, the improved PI controller designed by RA has better performance than Bode plot method.

SIMULATION AND EXPERIMENTAL RESULTS
To validate RA results, all the former case studies are simulated in MATLAB/Simulink. As discussed in Section 4, the inputs i d,re f and i q,re f can take any value in a circle centred at the origin with radius I nom . To emulate this condition, the PV unit is turned on and off with a switch every 0.4 s and its current i PV is either 0 or For experiments, the controllers are implemented using realtime simulator and controller OP5031 OPAL-RT. There is an I/O unit between OPAL-RT and VSC; it receives PWM signals from OPAL-RT and sends them to SKYPER 32 R gate drivers. It also receives the voltage and current measurements and sends them to OPAL-RT to be used in controller. Figure 8 shows the experimental setup. To provide random inputs for the experiments, the PV unit is connected for half of each case study time and the DC electronic load is randomly changed manually such that the d -and q-axes reference currents change within their limits. To implement the inductance error, three different inductors are used (0.5L f ,nom , L f ,nom , and 1.5L f ,nom ), while the controller always assumes L f ,nom .

System with inductance error
To evaluate transient behaviour of the grid-connected PV system, it is simulated for positive and negative step changes in dand q-axes reference currents and grid voltage. Figure 9 shows its transient response. The system operates in steady-state with i d,re f = 2 A, i q,re f = −2 A, and a randomly selected inductance as L f = 0.5L f,nom . Figure 9(  and IMC controller, respectively. The simulation results for the step changes in i d,re f and i q,re f show that the improved PI controller has a smaller rise time and settling time than the initial PI controller. Also, the IMC controller has the smallest rise time and settling time and it effectively decouples the d -and q-axes. Figure 9(c) shows the transient response when the grid voltage v s steps down to 0.95 pu at t = 0.5 s and steps up to 1 pu at t = 0.1 s. It shows that the system with all three controllers has robust performance during the grid voltage variation. The dand q-axes currents have lower oscillations in the system using improved PI and IMC controllers compared with the system using the initial PI controller. Figure 10 shows the simulation results for Cases A-I, A-II, and A-III discussed in Section 4. The maximum values of i t,d and i t,q are shown. Since the studied variable in this paper is ter- , RA results are presented as i t,d vs. i t,q . For consistency, simulation and experimental results are presented in similar trajectories. Figure 11 shows the simulation and experimental results for the d -and q-axis terminal currents. All results in each case are within the overapproximated reachable sets found for the respective cases in Section 4.1. Table 2 shows the maximum transient current in each case. As expected, the improved PI controller has a smaller transient current than the initial PI controller, and the IMC controller has the least transient current.
To compare the proposed improved PI controller and IMC controller with a robust control method, the studied PV system with the robust controller proposed in ref. [3] is simulated in MATLAB/Simulink. The designed robust controller parameters are q 1 = 10 16 , q 2 = 10 11.25 , q 4 = 10 10.5 , q 5 = 10 6 , and q 3 = q 6 = q 7 = q 8 = q 8 = 0. Figure 12 shows that the transient

FIGURE 10
Time-domain simulation results of the grid-connected PV unit with inductance error using (a) initial PI controller; (b) improved PI controller; (c) IMC controller current bound is 5.33 A for the system with robust controller. Considering Table 2, all three studied controllers have a transient current lower than the robust controller designed in [3]. Also, the simulated robust controller has nine tuned parameters which makes the design procedure more complicated than the improved PI controller (two parameters) and IMC controller (one parameter). Figure 13 shows the simulation and experimental results for the system without inductance error. The maximum transient current in each case study is shown in Table 2. The experimental current limits are higher than simulation results due to slight mismatch in parameters and nonidealities in measurements.

CONCLUSION
In this paper, the performance of a grid-connected PV unit is analysed using reachability analysis (RA). The reference terminal currents are uncertain-but-bounded inputs to the system