**3. Problem statement**

## **3.1. Network representation**

The present harmonic filtering study is based on the usual assumptions about the symmetry of the electrical network, the balance of the harmonic currents generated by the non-linear loads, and the independence between the harmonic orders. As a result, a single phase network is considered and the analysis is carried out for each individual harmonic order. Then, the relationship between nodal harmonic currents and voltages is defined through the admittance matrix [14], as follows:

$$\mathbf{I}\mathbf{\hat{}} = \mathbf{Y}\mathbf{\hat{}} \cdot \mathbf{V}\mathbf{\hat{}}\tag{3}$$

where

26 Simulated Annealing – Single and Multiple Objective Problems

**Figure 1.** SA algorithm

The driving mechanism of the SA process is described in Figure 1. The process of selecting and proposing a move is repeated until the system is considered in thermal equilibrium. Then, the temperature is reduced according to a temperature schedule and the system is allowed to reach thermal equilibrium again. Then, as the temperature of the system declines, the probability of accepting a worse move is decreased. This is the same as gradually moving to a frozen state in physical annealing. The process is finally stopped when no significant improvement is expected by further lowering temperature. At this point, the

The major advantage of a SA process over other methods is its ability to avoid becoming trapped in local minima. The algorithm employs a random search which not only accepts changes that decrease the objective function, but also some changes that can increase it. In consequence, a worse move can be accepted temporarily by leading subsequently to an improving solution. Then, the system sometimes goes uphill as well as downhill. However,

the lower the temperature, the less likely is any significant uphill excursion.

current state of the system is the solution to the optimisation problem.

*Ih* : the vector of the nodal harmonic currents,

*Yh* : the harmonic admittance matrix,

*Vh* : the vector of the nodal harmonic voltages,

*h* : the harmonic order.

The nodal admittance matrix is obtained from the impedance of every network component and the nodal currents are given from the harmonic currents generated by the non-linear loads.

In industrial distribution systems, loads are supplied through transformers by different voltage levels according to their rated power. Due to the voltage supplies, the filtering current is modified by the transformer ratio when the filter is connected to the primary or to the secondary side of the transformer. In consequence, the comparison between the set of solutions requires considering a per-unit system.

Let *Vn* be a vector composed of the nominal nodal voltages. The admittance matrix formulation (3) can be then rewritten in an equivalent system where the relationship (5) between nodal harmonic currents and voltages becomes independent from the nominal nodal voltages.

$$\mathbf{V}\_n \cdot \mathbf{I}\_h = \mathbf{V}\_n^2 \cdot \mathbf{Y}\_h \cdot \frac{\mathbf{V}\_h}{\mathbf{V}\_n} \tag{4}$$

$$\dot{\mathbf{i}}\_h = \mathbf{y}\_h \cdot \mathbf{v}\_h \tag{5}$$

Optimal Sizing of Harmonic Filters in Electrical Systems: Application of a Double Simulated Annealing Process 29

**Figure 2.** General flow chart of the optimisation procedure

specific application. It consists herein of the following input data:

*3.2.1. Initialisation procedure* 

The implementation of the algorithm is detailed further below and the focus is on the major parts of the procedure that can influence the behaviour of the SA approach to optimisation.

When thinking about any optimisation problem, one of the first considerations is to start with a suitable set of initial parameters in order to ensure a good algorithm performance, which means the quality of the solution returned and a reasonable computation time [15]. The initialisation procedure deals first with the definition of the problem's data linked to the

where

**i VI** *h nh* = ⋅ : the vector of the equivalent injected currents in kVA, *h nh* = ⋅ <sup>²</sup> **y VY** : the equivalent admittance matrix in kVA, *h h n* <sup>=</sup> **<sup>V</sup> v <sup>V</sup>** : the vector of the equivalent voltages in %.

Besides, an active filter connected to the node *k* can be considered as a current source which injects a current *jhk* on the network. A passive filter can be also modelled by a current source which magnitude would represent the harmonic current that must be drawn by the filter. In this last case, the phase angle depends on the harmonic voltage at the coupling node. All the additional currents produced by the set of filters are intended to reduce the harmonic voltages at each busbar. The location of these filters is defined in a node list called *ListFilter*. Assuming a given number of filters, the harmonic voltages are deduced from (6):

$$\mathbf{v}\_h = \mathbf{z}\_h \cdot (\mathbf{i}\_h \quad + \mathbf{j}\_h) \tag{6}$$

where

*jh* is the vector of the filtering current and *zh = yh-1* the equivalent impedance matrix.

#### **3.2. Optimisation process**

The aim of the current optimisation problem is to minimise the total filtering power to connect to the grid while the voltage standards are met for each nodal harmonic voltages. Prior to the search for the optimal current to be injected, it is necessary to answer a first question about the possibility to obtain a current vector *jh* able to reduce the voltages within the limiting levels *vlimit* at each node. For this reason, a double SA process is applied, as illustrated by the complete flow chart of figure 2. The purpose of the first one is to minimise the harmonic voltages, while checking if a filtering solution is existing or not.

In cases where no filtering solution can be found, the simulated voltage annealing is fully processed and returns the total filtering current *jtotal* able to check a minimum gap between the nodal harmonic voltages and their expected limits. Instead, when the voltage constraints are properly respected, the voltage annealing is partially processed. Once a local solution is achieved, the procedure is switched to a simulated current annealing which purpose is to minimise the total filtering current *jtotal* that meets the voltage requirements. The search space is actually the same regardless of the nature of the SA process; what is however different is the objective function, as will be argued further in the subsection 3.2.2.

**Figure 2.** General flow chart of the optimisation procedure

The implementation of the algorithm is detailed further below and the focus is on the major parts of the procedure that can influence the behaviour of the SA approach to optimisation.

## *3.2.1. Initialisation procedure*

28 Simulated Annealing – Single and Multiple Objective Problems

**i VI** *h nh* = ⋅ : the vector of the equivalent injected currents in kVA,

*h nh* = ⋅ <sup>²</sup> **y VY** : the equivalent admittance matrix in kVA,

**<sup>V</sup>** : the vector of the equivalent voltages in %.

where

*h*

where

**3.2. Optimisation process** 

*h*

*n* <sup>=</sup> **<sup>V</sup> v**

2 *h*

*n*

**<sup>V</sup>** ⋅= ⋅ ⋅ (4)

**i yv** *h hh* = ⋅ (5)

( ) **v zi j** *h hh h* =⋅ + (6)

*nh n h*

**<sup>V</sup> VI V Y**

Besides, an active filter connected to the node *k* can be considered as a current source which injects a current *jhk* on the network. A passive filter can be also modelled by a current source which magnitude would represent the harmonic current that must be drawn by the filter. In this last case, the phase angle depends on the harmonic voltage at the coupling node. All the additional currents produced by the set of filters are intended to reduce the harmonic voltages at each busbar. The location of these filters is defined in a node list called *ListFilter*.

Assuming a given number of filters, the harmonic voltages are deduced from (6):

*jh* is the vector of the filtering current and *zh = yh-1* the equivalent impedance matrix.

the harmonic voltages, while checking if a filtering solution is existing or not.

different is the objective function, as will be argued further in the subsection 3.2.2.

The aim of the current optimisation problem is to minimise the total filtering power to connect to the grid while the voltage standards are met for each nodal harmonic voltages. Prior to the search for the optimal current to be injected, it is necessary to answer a first question about the possibility to obtain a current vector *jh* able to reduce the voltages within the limiting levels *vlimit* at each node. For this reason, a double SA process is applied, as illustrated by the complete flow chart of figure 2. The purpose of the first one is to minimise

In cases where no filtering solution can be found, the simulated voltage annealing is fully processed and returns the total filtering current *jtotal* able to check a minimum gap between the nodal harmonic voltages and their expected limits. Instead, when the voltage constraints are properly respected, the voltage annealing is partially processed. Once a local solution is achieved, the procedure is switched to a simulated current annealing which purpose is to minimise the total filtering current *jtotal* that meets the voltage requirements. The search space is actually the same regardless of the nature of the SA process; what is however

When thinking about any optimisation problem, one of the first considerations is to start with a suitable set of initial parameters in order to ensure a good algorithm performance, which means the quality of the solution returned and a reasonable computation time [15]. The initialisation procedure deals first with the definition of the problem's data linked to the specific application. It consists herein of the following input data:


The initialisation procedure also starts with the parameters settings for the SA process, i.e. the cooling schedule that consists of the initial temperature and the rules for lowering it as the search for the optimal solution progresses. Each of these specific features will be introduced and explained, regarding the current issue about the optimal placement and sizing of harmonic filters in electrical systems.

• Starting temperature

In a typical SA process, the initial temperature is set sufficiently high to allow a move to almost any neighbourhood state. However, if the temperature starts at a too high value, the search may move to any neighbour and thus transform it into a random search, at least in the early stages. At the moment, there is no known method for finding a suitable starting temperature for a whole range of problems. An idea is to use the information on the cost function difference between one neighbour and another one to calculate a correct initial temperature. Another method suggested in [13] is to rapidly heat the system until a certain proportion of worse solutions are accepted and then slow cooling can start.

In the current application where a double SA process is considered, two initial temperatures for the voltage and the current annealing respectively are set as defined in (7) and (8). The both relationships mean that a 60% tolerance is allowed as regards the acceptance of an unfavourable harmonic voltage equal to twice the required limits and a jump of 30% in the amplitude of the average current is accepted with a same 60% tolerance.

$$Tv = -\frac{2 \times \max(\mathbf{v}\_{\text{limit}})}{\ln(0.6)}\tag{7}$$

Optimal Sizing of Harmonic Filters in Electrical Systems: Application of a Double Simulated Annealing Process 31

The temperature decreases during the search according to a function known as the cooling (or annealing) schedule. The way in which the temperature declines is critical to the success of the algorithm. Theory states that the number of iterations to execute at each temperature

In the literature, several theoretical and empirical control schemes are suggested [16-19] and can be categorised into classes such as monotonic schedules, geometric schedules, quadratic schedules and adaptive cooling schedules. Actually, many attempts have been made to derive or suggest good annealing schedules. Several comparative studies on the large

In a conventional SA process, the way to decrement the temperature is a simple linear method. By declining the temperature constantly, it provides the search with a higher transition probability in the beginning of the search and lower probability towards the end

The proposed annealing procedure tested on the minimum power sizing scheduling problem involves the above cooling schedule known as an exponential schedule (9), with a

The final decision to make is the number of iterations to consider at each temperature. A constant number seems to be an obvious scheme. An alternative is to dynamically change this number as the algorithm progresses. At lower temperature, it is advisable to set a large value so that the local optimum can be fully explored. At higher temperatures, the number

In the current issue, the number of steps *Nsteps* was set in connection to the number of

where *n* is a suitably large integer (*n* ≅ 200), defined to achieve the thermal equilibrium prior

equal to 0.95 in order to compromise between computation time and

should lie between [0.8 – 0.99], with better results found in the

α

( ) <sup>0</sup> *<sup>t</sup> Tt T* = ⋅*a* (9)

*N nN steps* =⋅ + ( 2) *filters* (10)

α (0 < α

, the longer it will take to

< 1). Then,

should be large enough to reach the thermal equilibrium.

variety of proposed cooling strategies are discussed in [20-22].

higher end of the range. However, the higher the value of

α

decrement the temperature to the stopping criterion.

of the search. An alternative is a geometric decrease by a constant factor

where *t* (for 'time') is the step count and *T0* is the starting temperature.

possible harmonic filters *Nfilters* to connect to the grid, as given in (10):

• Cooling schedule

experience shows that

α

of iterations can be less.

to a next temperature change.

optimisation performance.

• Iterations at each temperature

coefficient

$$T\_i = -\frac{0.3 \times \sum\_{q \in ListFilter} \left| \dot{h}^{hq} \right|}{\ln(0.6)}\tag{8}$$

• Final temperature

It is usual to let the temperature decrease until it reaches zero. However, this can make the algorithm run for too long, which should be a major drawback. In practice, the stopping criterion is a suitably low temperature, since the chances of accepting a worse move are almost the same when the temperature is null. In other words, the stopping criterion is met when the system is frozen at the current temperature; that is, no better or worse moves are being accepted.

Optimal Sizing of Harmonic Filters in Electrical Systems: Application of a Double Simulated Annealing Process 31

• Cooling schedule

30 Simulated Annealing – Single and Multiple Objective Problems

each node where a filter is connected,

sizing of harmonic filters in electrical systems.



The initialisation procedure also starts with the parameters settings for the SA process, i.e. the cooling schedule that consists of the initial temperature and the rules for lowering it as the search for the optimal solution progresses. Each of these specific features will be introduced and explained, regarding the current issue about the optimal placement and

In a typical SA process, the initial temperature is set sufficiently high to allow a move to almost any neighbourhood state. However, if the temperature starts at a too high value, the search may move to any neighbour and thus transform it into a random search, at least in the early stages. At the moment, there is no known method for finding a suitable starting temperature for a whole range of problems. An idea is to use the information on the cost function difference between one neighbour and another one to calculate a correct initial temperature. Another method suggested in [13] is to rapidly heat the system until a certain

In the current application where a double SA process is considered, two initial temperatures for the voltage and the current annealing respectively are set as defined in (7) and (8). The both relationships mean that a 60% tolerance is allowed as regards the acceptance of an unfavourable harmonic voltage equal to twice the required limits and a jump of 30% in the

> 2 max( ) limit ln(0.6)

<sup>0</sup> 0.3

*hq*

*j*

ln(0.6)

It is usual to let the temperature decrease until it reaches zero. However, this can make the algorithm run for too long, which should be a major drawback. In practice, the stopping criterion is a suitably low temperature, since the chances of accepting a worse move are almost the same when the temperature is null. In other words, the stopping criterion is met when the system is frozen at the current temperature; that is, no better or worse moves are

*q ListFilter*

å

*Tv* ´ = - **<sup>v</sup>** (7)

(8)


proportion of worse solutions are accepted and then slow cooling can start.

amplitude of the average current is accepted with a same 60% tolerance.

*i*

*T* <sup>Î</sup>

= -

´



• Starting temperature

• Final temperature

being accepted.

The temperature decreases during the search according to a function known as the cooling (or annealing) schedule. The way in which the temperature declines is critical to the success of the algorithm. Theory states that the number of iterations to execute at each temperature should be large enough to reach the thermal equilibrium.

In the literature, several theoretical and empirical control schemes are suggested [16-19] and can be categorised into classes such as monotonic schedules, geometric schedules, quadratic schedules and adaptive cooling schedules. Actually, many attempts have been made to derive or suggest good annealing schedules. Several comparative studies on the large variety of proposed cooling strategies are discussed in [20-22].

In a conventional SA process, the way to decrement the temperature is a simple linear method. By declining the temperature constantly, it provides the search with a higher transition probability in the beginning of the search and lower probability towards the end of the search. An alternative is a geometric decrease by a constant factor α (0 < α < 1). Then, experience shows that α should lie between [0.8 – 0.99], with better results found in the higher end of the range. However, the higher the value of α, the longer it will take to decrement the temperature to the stopping criterion.

The proposed annealing procedure tested on the minimum power sizing scheduling problem involves the above cooling schedule known as an exponential schedule (9), with a coefficient α equal to 0.95 in order to compromise between computation time and optimisation performance.

$$T(t) = T \circ \alpha^t \tag{9}$$

where *t* (for 'time') is the step count and *T0* is the starting temperature.

• Iterations at each temperature

The final decision to make is the number of iterations to consider at each temperature. A constant number seems to be an obvious scheme. An alternative is to dynamically change this number as the algorithm progresses. At lower temperature, it is advisable to set a large value so that the local optimum can be fully explored. At higher temperatures, the number of iterations can be less.

In the current issue, the number of steps *Nsteps* was set in connection to the number of possible harmonic filters *Nfilters* to connect to the grid, as given in (10):

$$N\_{\rm str\*} = n \cdot (N\_{\rm filters} + \mathcal{Z}) \tag{10}$$

where *n* is a suitably large integer (*n* ≅ 200), defined to achieve the thermal equilibrium prior to a next temperature change.
