# Discrete-Time Optimal Control of Photovoltaic Grid-Connected Inverter Based on Particle Swarm Optimization

Authors: Xinrui Liu; Guangru Zhang; Dongsheng Yang; Tongyu Shi; Xusheng He
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

This paper is concerned with the problem of optimal control of photovoltaic grid-connected inverter. Firstly, the discrete-time
nonlinear mathematical model of single-phase photovoltaic grid-connected inverter in the rotating coordinate system is constructed by the Delta operator, which simplifies the control process and facilitates direct digital realization. Then, a novel optimal control method which is significant to achieve trajectory
tracking for photovoltaic grid-connected inverter is developed by constructing a control Lyapunov functional against the difficulty of solving HJB (Hamilton-Jacobi-Bellman) partial differential equations. More importantly, the performance matrix of the controller which is very meaningful to many performance indicators of the system is also optimized via the particle swarm optimization. Finally, two numerical examples are provided to
illustrate the effectiveness of the proposed method.


---

## Body

## 1. Introduction

Energy is the basic material for survival and development of human, and it is also the blood of the world economy. Along with the development of society and economy and the improving of living standards, the demand for energy is increasingly urgent. In expected future, more than 60% of the worldwide energy will be indirectly consumed by conversion into electrical energy and thus it is urgent to take advantage of the new decentralized energy generation to meet the increasing demand for electricity [1, 2]. Considering many factors of the energy supply, solar energy is undoubtedly the ideal of sustainable development and green energy. Figure 1 shows a typical single-phase photovoltaic grid-connected system. The control of the system is composed by the boost chopper circuit and single-phase inverter. Boost chopper circuit is mainly to ensure photovoltaic modules running in MPPT (maximum power point tracking) point [3–7].Figure 1
Single-phase photovoltaic grid-connected system.In order to ensure control accuracy and obtain high-quality sine wave of the inverter output current which has the same frequency and phase with the grid voltage, experts and scholars have put forward many methods, such as the traditional PID control, dual-loop control, and hysteresis control. And repetitive control, proportional resonance control, and sliding mode control have also been proposed in recent years. Traditional methods [8–11] used to control the system are relatively simple and easy to adjust controller parameters, and their stability can be guaranteed and they have mature application, but it is difficult to adjust parameters of the system and eliminate static error and control accuracy is not high. In a word, the overall effect is not satisfactory. As the nonlinear characteristics of electronic devices and photovoltaic modules, it is difficult to obtain improved results by linear approximation using traditional methods to implement the control of the system, and these methods are more suitable for the system which can precisely obtain the steady-state point and always has small perturbations [12]. Therefore, it is essential to research nonlinear control methods for grid-connected photovoltaic inverter in order to obtain high-quality grid-connected inverter output current. In [13], it combines the parallel PI with repetitive control which solves the problem of phase lagging, but the improvement is limited and dynamic response can still be improved. In [14], using a Lyapunov functional to design the nonlinear controller for inverters has global stability, and its static and dynamic response get a better improvement, but its implementation is difficult and calculations are very complex. Sliding mode control having strong robustness and stability in a wide margin has been proved in [15, 16], but it cannot eliminate the steady-state error and also has oscillation problem. Optimal nonlinear control method using the power converter has been proposed in [17–19], and it has proved that the method can operate in a stable way in a wide fluctuating range and has strong robustness. Moreover, [17–20] also point out that the optimal control is more suitable for trajectory tracking problem, which also proves its stability and effectiveness, and [21, 22] also focus on the optimal control and prove that this method has good performance. Photovoltaic grid-connected inverter also can be considered in the current trajectory tracking problem, and using the optimal control can significantly improve performances of the system.However, it is worth noting that optimization problem generally solves the Hamilton-Jacobi-Bellman partial differential equations, which are usually difficult to solve. On the other hand, it is pointed out in [23, 24] that constructing a Lyapunov functional in the form of inverse optimality could overcome the shortcomings of optimal control. So far, due to the complexity of mathematical analysis, there has been less work on optimal control of photovoltaic grid-connected inverter. Therefore, this problem needs to be further investigated.In this paper, the problem of optimal control for photovoltaic grid-connected inverter is addressed. Firstly, the discrete-time nonlinear mathematical model of single-phase photovoltaic grid-connected inverter in the rotating coordinate system is constructed by the Delta operator in Section2, which simplifies the control process and facilitates direct digital realization. Then, a novel optimal control method for photovoltaic grid-connected inverter is given by constructing a control Lyapunov functional, and the optimization method for the performance matrix of the controller is also proposed in Section 3. Lastly, simulations are given to verify the above proposed method in Section 4.


## 2. Single-Phase Photovoltatic Grid-Connected Inverter Model

### 2.1. Transformation to Rotating Frame of Single Phase

Many parameters of power converters are AC variables which make the analysis of the system and the controller designer complex and cumbersome. Therefore, the introduction of the rotating coordinate system has an important value to simplify the variables.d-q rotating coordinate transformation is commonly used in the analysis and controller design of three-phase inverter. Shortcutting creating models of inverters and converters and all time-varying state variables of the topology have become constants, so that the analysis and design of the controller changes become even simpler [25]. Since d-q transformation requires at least two orthogonal phases, it is not useful for the single-phase inverter directly [26]. To construct a d-q rotating coordinate system of single-phase inverter, a virtual signal corresponding to the original signal is built to generate two orthogonal phase signals.Assume the actual steady-state variable is(1)XR=Xmsin(ωt+φ),
where Xm is the peak value of the variable, ω is the fundamental frequency, and φ is the initial phase of the sine variable. In ideal condition, the corresponding virtual orthogonal variable is
(2)XI=Xmcos⁡(ωt+φ).Transformation matrix from two-phase stationary coordinate to two-phase rotating coordinate can be obtained by using(3)T=[sin(ωt)cos⁡(ωt)cos⁡(ωt)-sin(ωt)].The variable in thed-q rotating coordinate system can be expressed as
(4)[XDXQ]=T[XRXI]=Xm[cos⁡φsinφ].Figure2 illustrates the entire building process from a stationary coordinate system to a rotating coordinate system. The time-varying variables of single-phase photovoltaic grid-connected inverter are mainly inverter output current and grid voltage. So the transformation of current and voltage are constructed in Figure 2, where uα represents the actual grid voltage and uβ represents the corresponding virtual orthogonal grid voltage; iα represents the actual inverter output current and iβ represents the corresponding virtual orthogonal inverter output current; uαβ represents the grid voltage and iαβ represents the inverter output current; ud represents the grid voltage in d axis and uq represents the grid voltage in q axis; id represents the inverter output current in d axis and iq represents the inverter output current in q axis. Figure 2 clearly shows that the time-varying variables of the inverter become constant variables, so that the analysis and design of the controller are simplified subsequently.Figure 2
Phasor diagram ofd-q coordinate.


### 2.2. The Mathematical Model of Single-Phase Grid-Connected Inverter

Figure1 shows a typical single-phase photovoltaic grid-connected inverter structure. From it, the mathematical model can be obtained as follows:
(5)U˙dc(t)=1CUdc(t)ugT(t)MPig(t),i˙g(t)=-AUdc(t)Lg-RgLgig(t)+2Udc(t)Lgu(t)-ug(t)Lg,
where ug=[ugdugq], ig=[igdigq], u=[uduq], A=[10], MP=[1001], and Udc is DC bus voltage, C is capacitor of DC bus, Lg is grid-connected filtering inductance, Rg is equivalent resistance of grid-connected filtering inductance, ugd and ugq are grid voltage in d axis and q axis, respectively, igd and igq are inverter output current in d axis and q axis, respectively, and ud and uq are the duty in d axis and q axis, respectively.Thanks to the low-cost microprocessor which is the main tool to implement the control method and requires the discrete samples for the signal processing normally, it is necessary to construct discrete-time model of the inverter which can realize the control strategy directly. Traditional shift operator is the main method to discretize continuous systems. When the sampling frequency of the discrete-time system increases, the traditional method is difficult to avoid the shortcomings of the following [27]. With the sampling frequency increasing, parameters of the discrete model are not inequivalent to the parameters of the corresponding continuous system, and when the relative degree of continuous system is greater than 1, some zero-polo will tend to the unit circle or go outside of the unit circle; on the other hand, the increasing sampling frequency will lead to the oscillation of the limit cycle and instability of the system. Therefore, when the sampling period is small, the traditional shift operator to realize discretization of systems will result in all poles locating on the boundary of stability for sampling systems and the stability will deteriotate. Control expert, Goodwin, proposes that Delta operator (incremental difference operator, often called δ operator) is a better method to discretize continuous systems [28, 29], and this method can make it close to the original continuous model in the fast sampling case and maintains the original kinetic properties. In this paper, δ operator is used to achieve the discrete-time mathematical model of the system.δ operator (incremental difference operator) is defined as
(6)δ(T)={dX(t)dt,(T=0),X(t+T)-X(t)T,(T≠0),
where T is sampling period of the system. When T=0, the system is continuous; when T≠0, the system is discrete.In order to design a discrete controller for digital implementation directly, the state space model can be obtained by theδ operator as follows:
(7)Udc(k+1)=Udc(k)+τ[1CUdc(k)ugT(k)MPig(k)],ig(k+1)=ig(k)+τ[-AUdc(k)Lg-RgLgig(k)ccccccccccc+2Udc(k)Lgu(k)-ug(k)Lg],
where τ is period of sample time. Obviously, the system matrix is with a nonlinear term. In Section 3, optimal control will be developed to deal with the output tracking problem of a general nonlinear system including (7).


## 2.1. Transformation to Rotating Frame of Single Phase

Many parameters of power converters are AC variables which make the analysis of the system and the controller designer complex and cumbersome. Therefore, the introduction of the rotating coordinate system has an important value to simplify the variables.d-q rotating coordinate transformation is commonly used in the analysis and controller design of three-phase inverter. Shortcutting creating models of inverters and converters and all time-varying state variables of the topology have become constants, so that the analysis and design of the controller changes become even simpler [25]. Since d-q transformation requires at least two orthogonal phases, it is not useful for the single-phase inverter directly [26]. To construct a d-q rotating coordinate system of single-phase inverter, a virtual signal corresponding to the original signal is built to generate two orthogonal phase signals.Assume the actual steady-state variable is(1)XR=Xmsin(ωt+φ),
where Xm is the peak value of the variable, ω is the fundamental frequency, and φ is the initial phase of the sine variable. In ideal condition, the corresponding virtual orthogonal variable is
(2)XI=Xmcos⁡(ωt+φ).Transformation matrix from two-phase stationary coordinate to two-phase rotating coordinate can be obtained by using(3)T=[sin(ωt)cos⁡(ωt)cos⁡(ωt)-sin(ωt)].The variable in thed-q rotating coordinate system can be expressed as
(4)[XDXQ]=T[XRXI]=Xm[cos⁡φsinφ].Figure2 illustrates the entire building process from a stationary coordinate system to a rotating coordinate system. The time-varying variables of single-phase photovoltaic grid-connected inverter are mainly inverter output current and grid voltage. So the transformation of current and voltage are constructed in Figure 2, where uα represents the actual grid voltage and uβ represents the corresponding virtual orthogonal grid voltage; iα represents the actual inverter output current and iβ represents the corresponding virtual orthogonal inverter output current; uαβ represents the grid voltage and iαβ represents the inverter output current; ud represents the grid voltage in d axis and uq represents the grid voltage in q axis; id represents the inverter output current in d axis and iq represents the inverter output current in q axis. Figure 2 clearly shows that the time-varying variables of the inverter become constant variables, so that the analysis and design of the controller are simplified subsequently.Figure 2
Phasor diagram ofd-q coordinate.


## 2.2. The Mathematical Model of Single-Phase Grid-Connected Inverter

Figure1 shows a typical single-phase photovoltaic grid-connected inverter structure. From it, the mathematical model can be obtained as follows:
(5)U˙dc(t)=1CUdc(t)ugT(t)MPig(t),i˙g(t)=-AUdc(t)Lg-RgLgig(t)+2Udc(t)Lgu(t)-ug(t)Lg,
where ug=[ugdugq], ig=[igdigq], u=[uduq], A=[10], MP=[1001], and Udc is DC bus voltage, C is capacitor of DC bus, Lg is grid-connected filtering inductance, Rg is equivalent resistance of grid-connected filtering inductance, ugd and ugq are grid voltage in d axis and q axis, respectively, igd and igq are inverter output current in d axis and q axis, respectively, and ud and uq are the duty in d axis and q axis, respectively.Thanks to the low-cost microprocessor which is the main tool to implement the control method and requires the discrete samples for the signal processing normally, it is necessary to construct discrete-time model of the inverter which can realize the control strategy directly. Traditional shift operator is the main method to discretize continuous systems. When the sampling frequency of the discrete-time system increases, the traditional method is difficult to avoid the shortcomings of the following [27]. With the sampling frequency increasing, parameters of the discrete model are not inequivalent to the parameters of the corresponding continuous system, and when the relative degree of continuous system is greater than 1, some zero-polo will tend to the unit circle or go outside of the unit circle; on the other hand, the increasing sampling frequency will lead to the oscillation of the limit cycle and instability of the system. Therefore, when the sampling period is small, the traditional shift operator to realize discretization of systems will result in all poles locating on the boundary of stability for sampling systems and the stability will deteriotate. Control expert, Goodwin, proposes that Delta operator (incremental difference operator, often called δ operator) is a better method to discretize continuous systems [28, 29], and this method can make it close to the original continuous model in the fast sampling case and maintains the original kinetic properties. In this paper, δ operator is used to achieve the discrete-time mathematical model of the system.δ operator (incremental difference operator) is defined as
(6)δ(T)={dX(t)dt,(T=0),X(t+T)-X(t)T,(T≠0),
where T is sampling period of the system. When T=0, the system is continuous; when T≠0, the system is discrete.In order to design a discrete controller for digital implementation directly, the state space model can be obtained by theδ operator as follows:
(7)Udc(k+1)=Udc(k)+τ[1CUdc(k)ugT(k)MPig(k)],ig(k+1)=ig(k)+τ[-AUdc(k)Lg-RgLgig(k)ccccccccccc+2Udc(k)Lgu(k)-ug(k)Lg],
where τ is period of sample time. Obviously, the system matrix is with a nonlinear term. In Section 3, optimal control will be developed to deal with the output tracking problem of a general nonlinear system including (7).


## 3. Controller Design

### 3.1. Mathematical Preliminaries of Optimal Control

From (7), the system can be considered as a nonlinear discrete-time affine system naturally. Therefore, optimal control commonly used in the discrete-time affine system can be used in the system. To understand the optimal control better, a brief discussion will be given in this section.Now, consider a typical nonlinear discrete-time affine system [30–33] as follows:
(8)xk+1=f(xk)+g(xk)uk,yk=h(xk),ccccccccccccccccccccccccccccccx0=x(0).
Along with the associated meaningful cost functional [34],
(9)V(xk)=∑n=k∞‍l(xn)+unTR(xn)un,
where xk∈Rn is the state of the system at time k∈N+={0,1,2,…}; uk∈Rm is the input of the control; f(xk):Rn→Rn, g(xk):Rn→Rn×m, and h(xk):Rn→Rn are smooth mappings; f(0)=0 and g(xk)≠0 for all xk≠0; V:Rn→R+; l:Rn→R+ is a positive semidefinite function (a function l(x) is a positive semidefinite function if, for all vectors, l(x)≥0; namely, there exist some vectors x for which l(x)=0 and for all others x,l(x)>0  [35]) and R:Rn→Rm×m is a real symmetric positive definite weighting matrix. In the section, we assume that all the states are measurable and available.And (9) can be transformed into
(10)V(xk)=l(xk)+ukTR(xk)uk+∑n=k+1∞‍(l(xn)+unTR(xn)un)=l(xk)+ukTR(xk)uk+V(xk+1),
where V(0)=0 is required as the boundary condition. Therefore, V(xk) becomes a Lyapunov functional.Based on Bellman’s optimality principle, the cost functionalV(xk) becomes invariant and satisfies the discrete-time Hamiltonian for the infinite horizon optimization case [34, 36], and it becomes
(11)V(xk)=min⁡{l(xk)+ukTR(xk)uk+V(xk+1)}.In order to obtain the control law of the system, (11) is transformed into
(12)Z(xk,uk)=l(xk)+ukTR(xk)uk+V(xk+1)-V(xk).The necessary condition to the optimal control law is(13)∂Z(xk,uk)∂uk=0.Combing (12) and (13), the determination of optimal control is corresponding to calculating the gradient of the right side of (12); then
(14)2R(xk)uk+∂V(xk+1)∂uk=0.Taking the derivative of (8) and substituting it into (14) yield
(15)2R(xk)uk+gT(xk)∂V(xk+1)∂xk+1=0.Then the optimal control law can be obtained by(16)-2R(xk)uk=gT(xk)∂V(xk+1)∂xk+1,
which can be rewritten as
(17)uk=-12R-1(xk)gT(xk)∂V(xk+1)∂xk+1.Substituting (17) into (11) gives
(18)V(xk)=l(xk)+14∂VT(xk+1)∂xk+1g(xk)[R-1(xk)]TR(xk)V(xk)=×R-1(xk)gT(xk)∂V(xk+1)∂xk+1+V(xk+1).Due to (ifA=AT, then (A-1)T=(AT)-1=A-1) [R-1(xk)]T=R-1(xk), one gets
(19)V(xk)=l(xk)+14∂VT(xk+1)∂xk+1g(xk)R-1(xk)V(xk)=×gT(xk)·∂V(xk+1)∂xk+1+V(xk+1),
which can be expressed as
(20)l(xk)+V(xk+1)-V(xk)+14∂VT(xk+1)∂xk+1g(xk)×R-1(xk)gT(xk)∂V(xk+1)∂xk+1=0.It is not simple to solve the HJB partial differential formulation (20) for V(xk). To overcome this awkward problem, inverse optimal control is proposed [37].The discrete-time inverse optimal control problem can be established as follows.Lemma 1 (see [37]).
Define the control law of the inverse optimal control(21)uk*=-12R-1(xk)gT(xk)∂V(xk+1)∂xk+1,
which is globally stabilizing if the following two points are satisfied:(1)
it achieves (global) stability whenxk=0 for system (8);(2)
V(xk) is a positive definite function (radial unbounded) and the following formulation is established:
(22)V-∶=V(xk+1)-V(xk)+uk*TR(xk)uk*≤0.
Since inverse optimal control is always used to get the minimum or maximumV(xk), a twice differentiable positive definite function can be proposed under meeting the requirements of Lemma 3.1 in [23]. Then,
(23)V(xk)=12xkTKTP-Kxk,
where P=KTP-K, P is a positive definite symmetric matrix, and K is an additional gain matrix to modify the convergence rate of the tracking error.
Substituting (23) into (21), it follows that
(24)uk*=-12R-1(xk)gT(xk)·[12×2KTP-Kxk+1]=-12R-1(xk)gT(xk)KTP-K[f(xk)+g(xk)uk*]=-12R-1(xk)gT(xk)KTP-Kf(xk)-12R-1(xk)gT(xk)KTP-Kg(xk)uk*.
Thenuk* can be formulated as
(25)uk*=-12[R(xk)+12gT(xk)KTP-Kg(xk)]-1×gT(xk)KTP-Kf(xk).


### 3.2. System Analysis

The control objective of single-phase grid-connected inverter generally is the output current. Now, it is hoped to feed energy to the grid with unity power factor in most cases, and in smart microgrid the inverter is required to compensate the reactive power according to the load conditions. Based on the above discussion, it is easy to know that the variables to be controlled are the DC bus voltage and reactive power feeding to the grid.The reactive power feeding to the grid can be defined as follows [38]:
(26)Qfg(k)=ugT(k)ig(k)2-(Upv(k)Ipv(k)ηug(k))2,
where ug(k) is the grid voltage, ig(k) is the inverter output current, Upv(k) and Ipv(k) are the voltage and current of the photovoltaic models, and η is the efficiency of the grid-connected device which is set to 1 in order to compute briefly.The static stable value using the reference for the reactive power can be expressed as(27)Qfgss(k)=Pfgss(k)fss1-fss2,
where Pfgss(k) is active power feeding to the grid and fss is static stable power factor.Therefore, the tracking errors of the state variables are defined as(28)eUdc(k)=Udc(k)-Udcss(k),eQfg(k)=Qfg(k)-Qfgss(k),
where Udc is DC bus voltage in real-time state, Udcss is DC bus voltage in static-steady state, eUdc is difference of the DC bus voltage, Qfg is reactive power feeding to the grid in real-time state, Qfgss is reactive power feeding to the grid in static-steady state, and eQfg is difference of the reactive power.Combining the discrete-time state space model of single-phase photovoltaic grid-connected inverter with the above system analysis, the steady-state operating point for reference can be formulated as(29)Udcss(k)+τCUdcss(k)ugT(k)MPigss(k)-Udcss(k+1)=0,ugT(k)igss(k)2-(Upv(k)Ipv(k)ug(k))2-ugT(k)MPigss(k)fss1-fss2=0.The current feeding to the grid and the DC bus voltage of steady-state operation point can be obtained by (29).


### 3.3. Optimal Controller

In order to simplify the controller design, the discrete-time state space model of single-phase photovoltaic grid-connected inverter should be expressed as a general nonlinear discrete-time affine system. Therefore, we define a synthesis vector composed by the DC bus voltage and the current feeding to the grid which is written asxk=[Udc(k)igd(k)igq(k)]T, and its corresponding steady-state vector is written as xkss=[Udcss(k)igdss(k)igqss(k)]T. Now, the model can be expressed as
(30)xk+1=f(xk)+g(xk)uk+hk,
where
(31)f(xk)=[Udc(k)+τCUdc(k)ugT(k)MPig(k)-τALgUdc(k)+(1-τRgLg)ig(k)],g(xk)=[00-2τLgUdc(k)00-2τLgUdc(k)],hk=[0-τugd(k)Lg-τugq(k)Lg].The tracking error of the state variables can be expressed as(32)exk=xk-xkss.Similarly,(33)exk+1=xk+1-xk+1ss.Combining (30), (32), and (33), we can get
(34)exk+1=f(exk)+f(xkss)+g(exk+xkss)uk+hk-xk+1ss.Generally, the composition of control variables is steady-state part and dynamic part. Therefore, we assumeuk=ukss+uk*, and (34) can be expressed as
(35)exk+1=f(exk)+f(xkss)+g(exk+xkss)ukss+g(exk+xkss)uk*+hk-xk+1ss.In order to make the formula match with the nonlinear discrete-time affine system, we assume that(36)f(xkss)+g(exk+xkss)ukss+hk-xk+1ss=0,g*(exk)=g(exk+xkss).Combining (35) and (36) yields
(37)ukss=(g(exk)+xkss)-1(-f(xkss)-hk+xk+1ss),(38)exk+1=f(exk)+g*(exk)uk*.Now, the steady-state part of the control variable can be obtained by (37). In addition, (38) is a typical nonlinear discrete-time affine system, so we can use the conclusion of Section 3.1. The dynamic part of the control variable can be formulated as
(39)uk*=-12[R(exk)+12g*T(exk)KTP-Kg*(exk)]-1uk*=g*T(exk)KTP-Kf(exk).The controller design is completed at this point.


### 3.4. Optimization of the Performance Matrix

Many performance indicators have relation with the choice of positive definite weighting matrix in the inverse optimal control; therefore it is very meaningful and important to optimize the definite weighting matrix of the optimal controller. In this paper, particle swarm algorithm is used to optimize the definite weighting matrix.Particle swarm optimization (PSO) is one of the intelligent optimization algorithms. In particle swarm, each particle can move at different speeds in the search space, and it is able to remember the optimal position which has been visited. In addition, all groups can share all optimal locations of the individual particles. Each particle and group can mediate the movement speed and direction according to a certain strategy and their best position, and the search space tends to be the optimal position for solving ultimately [39]. Each particle can be updated according to the following formula in the solution space of their own position and velocity:
(40)Vidk+1=w(k)Vidk+c1r1(Pidk-Xidk)+c2r2(Pgdk-Xidk),Xidk+1=Xidk+Vidk+1,w(k)=wstart(wstart-wend)(Tmax⁡-k)Tmax⁡,
where w(k) is the inertia weight; wstart is the initial inertia weight; wend is the inertia weight of the iteration to the maximum number of times; k is the number for the current iteration; Tmax⁡ is the maximum number of iterations; d=1,2,…,D; k=1,2,…,n; the entire group has n particles, and each particle is a vector with D-dimension; Xid is the position of particle; Vid is the velocity of particle; Pid is extreme value of individual; Pgd is extreme value for the group; c1 and c2 are nonnegative constants, called acceleration factor; r1 and r2 are random intervals number distributed in [0,1].By optimizing, the fitness function for particle is defined as follows:(41)fpso=∑i=1M‍((∑k=1N‍|Δx(i,k)|)/N)M,
in which M is number of state variables to be controlled which is equal to the dimension of the search space; N is the size of population.


## 3.1. Mathematical Preliminaries of Optimal Control

From (7), the system can be considered as a nonlinear discrete-time affine system naturally. Therefore, optimal control commonly used in the discrete-time affine system can be used in the system. To understand the optimal control better, a brief discussion will be given in this section.Now, consider a typical nonlinear discrete-time affine system [30–33] as follows:
(8)xk+1=f(xk)+g(xk)uk,yk=h(xk),ccccccccccccccccccccccccccccccx0=x(0).
Along with the associated meaningful cost functional [34],
(9)V(xk)=∑n=k∞‍l(xn)+unTR(xn)un,
where xk∈Rn is the state of the system at time k∈N+={0,1,2,…}; uk∈Rm is the input of the control; f(xk):Rn→Rn, g(xk):Rn→Rn×m, and h(xk):Rn→Rn are smooth mappings; f(0)=0 and g(xk)≠0 for all xk≠0; V:Rn→R+; l:Rn→R+ is a positive semidefinite function (a function l(x) is a positive semidefinite function if, for all vectors, l(x)≥0; namely, there exist some vectors x for which l(x)=0 and for all others x,l(x)>0  [35]) and R:Rn→Rm×m is a real symmetric positive definite weighting matrix. In the section, we assume that all the states are measurable and available.And (9) can be transformed into
(10)V(xk)=l(xk)+ukTR(xk)uk+∑n=k+1∞‍(l(xn)+unTR(xn)un)=l(xk)+ukTR(xk)uk+V(xk+1),
where V(0)=0 is required as the boundary condition. Therefore, V(xk) becomes a Lyapunov functional.Based on Bellman’s optimality principle, the cost functionalV(xk) becomes invariant and satisfies the discrete-time Hamiltonian for the infinite horizon optimization case [34, 36], and it becomes
(11)V(xk)=min⁡{l(xk)+ukTR(xk)uk+V(xk+1)}.In order to obtain the control law of the system, (11) is transformed into
(12)Z(xk,uk)=l(xk)+ukTR(xk)uk+V(xk+1)-V(xk).The necessary condition to the optimal control law is(13)∂Z(xk,uk)∂uk=0.Combing (12) and (13), the determination of optimal control is corresponding to calculating the gradient of the right side of (12); then
(14)2R(xk)uk+∂V(xk+1)∂uk=0.Taking the derivative of (8) and substituting it into (14) yield
(15)2R(xk)uk+gT(xk)∂V(xk+1)∂xk+1=0.Then the optimal control law can be obtained by(16)-2R(xk)uk=gT(xk)∂V(xk+1)∂xk+1,
which can be rewritten as
(17)uk=-12R-1(xk)gT(xk)∂V(xk+1)∂xk+1.Substituting (17) into (11) gives
(18)V(xk)=l(xk)+14∂VT(xk+1)∂xk+1g(xk)[R-1(xk)]TR(xk)V(xk)=×R-1(xk)gT(xk)∂V(xk+1)∂xk+1+V(xk+1).Due to (ifA=AT, then (A-1)T=(AT)-1=A-1) [R-1(xk)]T=R-1(xk), one gets
(19)V(xk)=l(xk)+14∂VT(xk+1)∂xk+1g(xk)R-1(xk)V(xk)=×gT(xk)·∂V(xk+1)∂xk+1+V(xk+1),
which can be expressed as
(20)l(xk)+V(xk+1)-V(xk)+14∂VT(xk+1)∂xk+1g(xk)×R-1(xk)gT(xk)∂V(xk+1)∂xk+1=0.It is not simple to solve the HJB partial differential formulation (20) for V(xk). To overcome this awkward problem, inverse optimal control is proposed [37].The discrete-time inverse optimal control problem can be established as follows.Lemma 1 (see [37]).
Define the control law of the inverse optimal control(21)uk*=-12R-1(xk)gT(xk)∂V(xk+1)∂xk+1,
which is globally stabilizing if the following two points are satisfied:(1)
it achieves (global) stability whenxk=0 for system (8);(2)
V(xk) is a positive definite function (radial unbounded) and the following formulation is established:
(22)V-∶=V(xk+1)-V(xk)+uk*TR(xk)uk*≤0.
Since inverse optimal control is always used to get the minimum or maximumV(xk), a twice differentiable positive definite function can be proposed under meeting the requirements of Lemma 3.1 in [23]. Then,
(23)V(xk)=12xkTKTP-Kxk,
where P=KTP-K, P is a positive definite symmetric matrix, and K is an additional gain matrix to modify the convergence rate of the tracking error.
Substituting (23) into (21), it follows that
(24)uk*=-12R-1(xk)gT(xk)·[12×2KTP-Kxk+1]=-12R-1(xk)gT(xk)KTP-K[f(xk)+g(xk)uk*]=-12R-1(xk)gT(xk)KTP-Kf(xk)-12R-1(xk)gT(xk)KTP-Kg(xk)uk*.
Thenuk* can be formulated as
(25)uk*=-12[R(xk)+12gT(xk)KTP-Kg(xk)]-1×gT(xk)KTP-Kf(xk).


## 3.2. System Analysis

The control objective of single-phase grid-connected inverter generally is the output current. Now, it is hoped to feed energy to the grid with unity power factor in most cases, and in smart microgrid the inverter is required to compensate the reactive power according to the load conditions. Based on the above discussion, it is easy to know that the variables to be controlled are the DC bus voltage and reactive power feeding to the grid.The reactive power feeding to the grid can be defined as follows [38]:
(26)Qfg(k)=ugT(k)ig(k)2-(Upv(k)Ipv(k)ηug(k))2,
where ug(k) is the grid voltage, ig(k) is the inverter output current, Upv(k) and Ipv(k) are the voltage and current of the photovoltaic models, and η is the efficiency of the grid-connected device which is set to 1 in order to compute briefly.The static stable value using the reference for the reactive power can be expressed as(27)Qfgss(k)=Pfgss(k)fss1-fss2,
where Pfgss(k) is active power feeding to the grid and fss is static stable power factor.Therefore, the tracking errors of the state variables are defined as(28)eUdc(k)=Udc(k)-Udcss(k),eQfg(k)=Qfg(k)-Qfgss(k),
where Udc is DC bus voltage in real-time state, Udcss is DC bus voltage in static-steady state, eUdc is difference of the DC bus voltage, Qfg is reactive power feeding to the grid in real-time state, Qfgss is reactive power feeding to the grid in static-steady state, and eQfg is difference of the reactive power.Combining the discrete-time state space model of single-phase photovoltaic grid-connected inverter with the above system analysis, the steady-state operating point for reference can be formulated as(29)Udcss(k)+τCUdcss(k)ugT(k)MPigss(k)-Udcss(k+1)=0,ugT(k)igss(k)2-(Upv(k)Ipv(k)ug(k))2-ugT(k)MPigss(k)fss1-fss2=0.The current feeding to the grid and the DC bus voltage of steady-state operation point can be obtained by (29).


## 3.3. Optimal Controller

In order to simplify the controller design, the discrete-time state space model of single-phase photovoltaic grid-connected inverter should be expressed as a general nonlinear discrete-time affine system. Therefore, we define a synthesis vector composed by the DC bus voltage and the current feeding to the grid which is written asxk=[Udc(k)igd(k)igq(k)]T, and its corresponding steady-state vector is written as xkss=[Udcss(k)igdss(k)igqss(k)]T. Now, the model can be expressed as
(30)xk+1=f(xk)+g(xk)uk+hk,
where
(31)f(xk)=[Udc(k)+τCUdc(k)ugT(k)MPig(k)-τALgUdc(k)+(1-τRgLg)ig(k)],g(xk)=[00-2τLgUdc(k)00-2τLgUdc(k)],hk=[0-τugd(k)Lg-τugq(k)Lg].The tracking error of the state variables can be expressed as(32)exk=xk-xkss.Similarly,(33)exk+1=xk+1-xk+1ss.Combining (30), (32), and (33), we can get
(34)exk+1=f(exk)+f(xkss)+g(exk+xkss)uk+hk-xk+1ss.Generally, the composition of control variables is steady-state part and dynamic part. Therefore, we assumeuk=ukss+uk*, and (34) can be expressed as
(35)exk+1=f(exk)+f(xkss)+g(exk+xkss)ukss+g(exk+xkss)uk*+hk-xk+1ss.In order to make the formula match with the nonlinear discrete-time affine system, we assume that(36)f(xkss)+g(exk+xkss)ukss+hk-xk+1ss=0,g*(exk)=g(exk+xkss).Combining (35) and (36) yields
(37)ukss=(g(exk)+xkss)-1(-f(xkss)-hk+xk+1ss),(38)exk+1=f(exk)+g*(exk)uk*.Now, the steady-state part of the control variable can be obtained by (37). In addition, (38) is a typical nonlinear discrete-time affine system, so we can use the conclusion of Section 3.1. The dynamic part of the control variable can be formulated as
(39)uk*=-12[R(exk)+12g*T(exk)KTP-Kg*(exk)]-1uk*=g*T(exk)KTP-Kf(exk).The controller design is completed at this point.


## 3.4. Optimization of the Performance Matrix

Many performance indicators have relation with the choice of positive definite weighting matrix in the inverse optimal control; therefore it is very meaningful and important to optimize the definite weighting matrix of the optimal controller. In this paper, particle swarm algorithm is used to optimize the definite weighting matrix.Particle swarm optimization (PSO) is one of the intelligent optimization algorithms. In particle swarm, each particle can move at different speeds in the search space, and it is able to remember the optimal position which has been visited. In addition, all groups can share all optimal locations of the individual particles. Each particle and group can mediate the movement speed and direction according to a certain strategy and their best position, and the search space tends to be the optimal position for solving ultimately [39]. Each particle can be updated according to the following formula in the solution space of their own position and velocity:
(40)Vidk+1=w(k)Vidk+c1r1(Pidk-Xidk)+c2r2(Pgdk-Xidk),Xidk+1=Xidk+Vidk+1,w(k)=wstart(wstart-wend)(Tmax⁡-k)Tmax⁡,
where w(k) is the inertia weight; wstart is the initial inertia weight; wend is the inertia weight of the iteration to the maximum number of times; k is the number for the current iteration; Tmax⁡ is the maximum number of iterations; d=1,2,…,D; k=1,2,…,n; the entire group has n particles, and each particle is a vector with D-dimension; Xid is the position of particle; Vid is the velocity of particle; Pid is extreme value of individual; Pgd is extreme value for the group; c1 and c2 are nonnegative constants, called acceleration factor; r1 and r2 are random intervals number distributed in [0,1].By optimizing, the fitness function for particle is defined as follows:(41)fpso=∑i=1M‍((∑k=1N‍|Δx(i,k)|)/N)M,
in which M is number of state variables to be controlled which is equal to the dimension of the search space; N is the size of population.


## 4. Simulation Results

The performances of the proposed nonlinear optimal controller are illustrated by simulation using MATLAB. The principle of simulation setup based on the experimental schematic Figure1 is described by Figure 3.Figure 3
The simulation flow diagram.In order to verify that the controller has strong robustness and weak dependence on the parameters of the system, two simulations will be given. In the two simulations, the parameters of particle swarm are set as follows: iterationTmax⁡ is 1000, particle d is 100, c1 and c2 are 2, wstart is 0.9, and wend is 0.4.(1)  The First Simulation. The control design parameters and values of the first simulation are given in Table 1.Table 1
Characteristics of the controlled system.
Parameter
Symbol
Value
(1)
DC link capacitor
C
dc
4700μF
(2)
Grid filter inductor
L
g
2.1 mH
R
g
0.6 Ω
(3)
Switching frequency
τ
50 kHz
(4)
DC link voltage
U
dc
40 V
(5)
Grid voltage
u
g
22 V
(6)
Grid frequency
ω
50 HzThe objection function evaluation by PSO is shown in Figure4.Figure 4
Best individual fitness by performance optimization.Optimized performance matrix is(42)P=[0.75720.36850.61770.36850.73020.23620.61770.23620.7487].Dynamic tracking of the DC bus voltage response curve is shown in Figure5, wherein the red means the dynamic tracking state values and blue indicates its reference value. Since grid is equivalent to infinite load in the system, its d axis reference current is set to 0. In order to ensure the unity power factor, the q axis current reference value is also set to 0 apparently. Thus, the dynamic current tracking response curve of d axis and q axis is shown in Figures 6 and 7, similarly.Figure 5
The response curve of DC bus voltage.Figure 6
The response curve ofd axis current.Figure 7
The response curve ofq axis current.(2)  The Second Simulation. The control design parameters and values of the second simulation are given in Table 2.Table 2
Characteristics of the controlled system.
Parameter
Symbol
Value
(1)
DC link capacitor
C
dc
6800μF
(2)
Grid filter inductor
L
g
1 mH
R
g
0.32 Ω
(3)
Switching frequency
τ
45 kHz
(4)
DC link voltage
U
dc
40 V
(5)
Grid voltage
u
g
22 V
(6)
Grid frequency
ω
50 HzThe objection function evaluation by PSO is shown in Figure8.Figure 8
Best individual fitness by performance optimization.Optimized performance matrix can be obtained like simulation 1, so it is not necessary to repeat it here.Dynamic tracking of the DC bus voltage response curve is shown in Figure9, wherein the red means the dynamic tracking state values and blue indicates its reference value. Since grid is equivalent to infinite load, its d axis reference current is set to 0. In order to ensure the unity power factor, the q axis current reference value is also set to 0 apparently. Thus, the dynamic current tracking response curve of d axis and q axis is shown in Figures 10 and 11, similarly.Figure 9
The response curve of DC bus voltage.Figure 10
The response curve ofd axis current.Figure 11
The response curve ofq axis current.From Figures5, 6, 7, 9, 10, and 11, it is obvious that the inverse optimal controller can realize zero steady-state error and has good dynamic and steady response characteristics. Comparing the results of simulation 1 and simulation 2, it is easy to prove that the controller has strong robustness to the changing parameters of the system. In a word, the designed controller has less dependence on the system parameters and has good performance.


## 5. Conclusion

In order to achieve the active and reactive power independently and use the low-cost microcontroller to implement the control strategy directly, a discrete-time model of single-phase photovoltaic grid-connected inverter in the rotating coordinate system has been constructed in this paper. Since the bandwidth of the controller is required not to be very high and it has good tracking performance of the optimal control, this paper has proposed an optimal controller for the single-phase inverter. Against the difficulty of solving the HJB partial differential equations, a Lyapunov function has been constructed to solve the problem perfectly. In addition, the performance matrix has been optimized by the particle swarm optimization. The simulation results have shown that the proposed method is better to eliminate static error of the system, has good response characteristic, and has good robustness and following performance. The future researches are to study the optimal control for the whole systems with the preboost conversion portion and study the coordinate control of the active and reactive power.



---

*Source: 249312-2014-05-20.xml*