Solving Optimization Problems Using PSO algorithm.

Introduction

As we have discussed the difference between the conventional optimization algorithms and the unconventional algorithms are summarized below.

1.The conventional algorithms start from a single  initial feasible  solution where as the unconventional algorithms start with a set of initial(feasibility not a prerequisite)  solutions. 

2.The conventional algorithms need the function to be continuously differentiable through out the range of search.

3.The accuracy of the conventional algorithms depending on the selection of initial solution.Always a local solution is assured. The performance of the unconventional algorithms will  vary for every run. A solution is assured.It may be local or global.Very rarely diverges.

Particle Swarm Optimization(PSO)

This  PSO algorithm also one of the important unconventional optimization algorithms.PSO optimizes a problem by having a population of candidate solutions, here dubbed particles, and moving these particles around in the search-space according to simple mathematical formulae. The movements of the particles are guided by the best found positions in the search-space which are updated as better positions are found by the particles.

PSO is originally attributed to Kennedy, Eberhart and Shi  and was first intended for simulating social behaviour. The algorithm was simplified and it was observed to be performing optimization. The book by Kennedy and Eberhart  describes many philosophical aspects of PSO and swarm intelligence.

As I am more interested in the implementation of this algorithm interested readers can get more details in the web regrading this algorithm.To implement PSO you neeed not know anything about the algorithm,you should know how to use the code to solve your modelled problem.

http://en.wikipedia.org/wiki/Particle_swarm_optimization

PSO Code in MATLAB

There are so many variants of this PSO the code  which  I found simple and powerful is the vectorised PSO  Toolbox by Mr Brian Birge.I thank him for the excellent toolbox.It can be downloaded from the link.

http://www.mathworks.com/matlabcentral/fileexchange/7506

Just download it and unzip it as a folder.You can either add it to MATLAB path or make it as  the default folder.

Solving an Optimization problem by PSO

For implementing  most of  the unconventional algorithms(including the PSO)   follow the steps.

1. Model the problem as unconstrained minimization problem( you can refer my earlier post on MATLAB GA toolbox ).

2. Write it as a function file.There are two types of functiones (i) scalar function(ii)vector function

In the scalar function for for every value of X the corresponding function value is returned.

X1-----------F(X1)

In the vector function for for every set of X values the corresponding function values are returned.

[X1 X2 X3 X4...Xn]-------------[F(X1),F(X2),F(X3),F(X4)...F(Xn)]

The code by Mr Brian Birge requires a vector function file .

3.After writing the function file save it in the same PSO folder and change it as default.

The toolbox has lot of details regarding the setting..

4.Run it .You get the solution as well the progressive graphs( iteration vs best solution and movement of particles in two dimension) .

Implementation

The same problem discussed in GA tool box is considered. Same way two files have to be written the first one is  the pso setting file(test1.m) .This calls the subroutine vector function file ex2.m.

program 1 test1.m

clear;
clc
% lower limir
l=[0 0];
% upper limit
u=[ 10 10];
ran=[l' u'];
% number of variables
n=2;
% settings for pso
%1. population,2. no of iterations,3 refresh on screen 4&5 type of PSO 06&7
%are final and initial swarm velocity 8. final iteration to reach final
%velocity 9. expected function value
Pdef = [20 200 10 2 2 0.9 0.4 200 50 5000 NaN 0 0];
[OUT]=pso_Trelea_vectorized('ex2',n,1,ran,0,Pdef);
out=abs(OUT)
P=out(1:n)

[F x]=ex2(P')

Program 2 vector function file ex2.m

function [F x]=ex2(x);

F=sum((x.*x).').'+1000*((sum(x')-10))^2;

Results

PSO: 1/200 iterations, GBest =                  100.

PSO: 20/200 iterations, GBest =   52.210493585432118.

PSO: 40/200 iterations, GBest =   50.166913798388642.

PSO: 60/200 iterations, GBest =   49.997052264769565.

PSO: 80/200 iterations, GBest =   49.978589042793075.

PSO: 100/200 iterations, GBest =   49.978455656591841.

PSO: 120/200 iterations, GBest =   49.975851601760063.

PSO: 140/200 iterations, GBest =   49.975342456412399.

PSO: 160/200 iterations, GBest =   49.975036994685418.

PSO: 180/200 iterations, GBest =    49.97501323855063.

PSO: 200/200 iterations, GBest =   49.975012959884864.

out =

   4.99798359369474

   4.99701881234046

  49.97501295916326

P =

   4.99798359369474

   4.99701881234046

F =

  49.97501295916326

x =

   4.99798359369474   4.99701881234046

You are welcome to make some suggestions .All the best.

Power System State Estimation: Part 3

If G1 is having a rank of ‘n’ then the system is fully observable.

Algorithm

1. Start iterations, set the iteration index i=0.

2. Initialize the state vector x₀, typically as a flat start.

3. Calculate the matrix G(x₀) as per (4)

4. Calculate the matrix H as per (5)

5. Determine x₁ as per (6).

6. Test for convergence, max [G(x₁)]< e

7. If no, update x , and go to step 3. If yes stop.

For power system state estimation the matrix H is the jacobian matrix of the Power flow the difference is that the jacobian is not square it is mX n

Power System State Estimation: Part 2

Least Squares Estimation & Weighted Least Squares Estimation

If the number of equations(linear or non linear ) is more than the no of variables to be determined least squares estimation is used to find a compromising solution.First let us discuss the linear least square problem ,then we turn our attention to the non linear and weighted least square problems.

The linear least square problem is determined as follows.

Determine X
 which satisfy the set of equations

A*X=b
                                                                                                                                         (1)

where A is a m X n matrix where m is the number of linear equations.

X =[x1 x2 ...xn] the solution vector of nX1 ,where m>n

so A is not a square matrix it's rank may be equal or less than 'n'.  

To solve this system of equations the concept of least squares is used.

The variable J is the summation of least squares. The solution vector which minimize the  ‘’J is called least squares solution.

J=(A*X-b)^T*(A*X-b)                                                                                                                        (2)

By applying khun tucker conditions the condition for optimality can be derived as .

A^T*A*X=A^T*b                                                                                                                                  (3)

 If the term A^T*A is having rank of ‘n’ the solution is easily determined  bys solving (3).If it is less than ‘n’ Orthogonal-triangular decomposition. is used to solve the system of equations.





% Example 1. Start with
%program to solve linear least squares pproblem
clear;
clc;
A =  [ 1     2     3
       4     5     6
       7     8     9
      10    11    12 ]
% This is a rank-deficient matrix; the middle column is the average of the other two columns. The rank deficiency is revealed by the factorization: [Q,R] = qr(A)




b = [1;3;5;7]
[Q,R] = qr(A)
  x = R\(R'\(A'*b))
        r = b - A*x
        e = R\(R'\(A'*r))
        x = x + e;


Solution

A =


     1     2     3
     4     5     6
     7     8     9
    10    11    12




b =


     1
     3
     5
     7




Q =


   -0.0776   -0.8331    0.5456   -0.0478
   -0.3105   -0.4512   -0.6919    0.4704
   -0.5433   -0.0694   -0.2531   -0.7975
   -0.7762    0.3124    0.3994    0.3748




R =


  -12.8841  -14.5916  -16.2992
         0   -1.0413   -2.0826
         0         0   -0.0000
         0         0         0


Warning: Rank deficient, rank = 2,  tol =   2.2550e-014.


x =


    0.5000
         0
    0.1667




r =


  1.0e-013 *


    0.3575
    0.1732
   -0.0089
   -0.1954


Warning: Rank deficient, rank = 2,  tol =   2.2550e-014.
e =


  1.0e-013 *


   -0.2709
         0
    0.2095

Power System State Estimation: Part 1

Introduction

Power system state estimation is defined as the act of estimating the state of the network from the redundant telemetry measurements. Static state estimation refers to the procedure of obtaining the voltage phasors at all of the system buses at a given point in time. This can be achieved by direct means which involve very accurate synchronized phasor measurements of all bus voltages in the system. However, such an approach would be very vulnerable to measurement errors or telemetery failures. 

Instead,state estimation procedure makes use of a set of redundant measurements in order to filter out such errors and find an optimal estimate.The measurements may include not only the conventional power and voltage measurements, but also those others such as the current magnitude or synchronized voltage phasor measurements as well. Simultaneous measurement of quantities at different parts of the system is practically impossible, hence a certain amount of time skew between measurements is commonly tolerated. This tolerance is justified due to the slowly varying operating conditions of the power systems under normal operating conditions..

Data for State Estimation

1. Network data

2. Measured Data

(i) Real Power (P)

(ii) Reactive Power (Q)

(iii) Real Line Flow (P_ij)

(iv) Reactive Line Flow (Q_ij)
(v)Magnitude of Line current( [I_ij] )

(vi) Voltage magnitude ([V])

Solving Linear and Quadratic Programming Problems by MATLAB

Introduction

Optimization is defined as Minimizing (or Maximizing) an objective function subject to some constraints .If the objective function and the all the constrains are linear it is called linear programming. If the objective function is quadratic and the all the constraints are linear, it is known as quadratic programming. Let us learn how to use the linear and quadratic programming routines of the MATLAB.

Linear programming

The command for implementing the matlab linear programming routine is ‘linprog’.

Let us consider a linear programming problem with ‘n’ variables, ‘m’ inequality constraints, ‘k’ equality constraints, and the lower ,upper bounds LB and UB.

The linear programming is defined in matlab like this

min f'*x    subject to:   A*x <= b, Aeq*x = beq ,LB≤ x ≤ UB

x-=[x₁ x₂ x₃ ….x_n]^T x is a column matrix of 1 X n) to be determined.  

f=[c₁ c₂ c₃ ….c_n]  f is a vector(row matrix of n X 1)

A is the matrix of inequalities having a size of m X n

Aeq is the matrix of equalities having a size of k X n

LB-=[l₁ l₂ l₃ ….l_n]^T x is a column matrix of 1 X n)

UB= [u₁ u₂, u₃ ….u_n]^T x is a column matrix of 1 X n)

Example

f=[10 5 12 13 20];

A= [10     6     1     0     8

     2     8     4     7     0

     6     9     8     4     7

     5     7     0     9     4

     9     2     1     5     8

     8     4     2     4     5

     5     9     2     8     7

     0     9     6     5     4

     8     4     3     2     3

     4     9     2     7     2];

b=[70    63   104    80    81    68   101    81    53    71]';

Aeq=[1 1 1 1 1]

beq=15;

LB=[0 0 0 0 0]'

UB=[20 20 20 20 20]'

X=linprog(f,A,b,Aeq,beq,LB,UB)

Solution

Optimization terminated.

X =

    2.1683

    0.0000

    9.6931

    2.8416

    0.2970

Quadratic programming

The command for implementing the matlab quadratic programming routine is ‘quadprog’.Let us consider a quadratic programming problem with ‘n’ variables, ‘m’ inequality constraints, ‘k’ equality constraints, and the lower, upper bounds LB and UB.

The quadratic programming is defined in matlab like this

   min 0.5*x'*H*x + f'*x    subject to:   A*x <= b, Aeq*x = beq ,LB≤ x ≤ UB

x-=[x₁ x₂ x₃ ….x_n]^T x is a column matrix of 1 X n) to be determined.  

H= is a square matrix of n X n

f=[c₁ c₂ c₃ ….c_n]  f is a vector(row matrix of n X 1)

A is the matrix of inequalities having a size of m X n

Aeq is the matrix of equalities having a size of k X n

LB-=[l₁ l₂ l₃ ….l_n]^T x is a column matrix of 1 X n)

UB= [u₁ u₂, u₃ ….u_n]^T x is a column matrix of 1 X n)

Example

Minimize x₁^2+x₁x₂+x₂₂+3x₁+5x₂

Subject to the constraints

x₁+x₂=20;

2x₁+3x₂≤100

[0 0]^T≤x₁,x₂≤[50 50]^T

H=[2 1;1 2];

f=[3 5];

A= [2 3];

b=[100];

Aeq=[1 1]

beq=20;

LB=[0 0 0 0 0]';

UB=[50 50]';

[X ff]=quadprog(H,f,A,b,Aeq,beq,LB,UB)

Solution

X =

   11.0000

    9.0000

ff is the objective function value

ff =  379.0000

Power Flow Analysis: Transmission and Distribution Systems: An Introduction

Power flow or load flow is one of the important power system analyses. Problem. This problem is determining the unspecified (dependent) parameters of the power network by some specified (independent) parameters.

It is known that the power system consists of generators, loads, transformers, transmission lines and other protective devices. In power flow every bus (node) is having four parameters.

1. Real Power P

2. Reactive power Q

3. Voltage magnitude [V]

4. Voltage magnitude φ

Depending on the type of the bus two out the four variables are specified and the remaining two has to be determined. The specified two variables and unspecified two variables of each bus are related by the kirchoff’s network laws.

The network elements are modeled according to their steady state characteristics.

1. Generator –voltage source with series impedance

2. Loads: constant power or constant admittance (impedance)

3. Transmission line: pi model

4. Transformers are to be incorporated in the transmission line model.

The specified and unspecified parameters
Type of the bus	Specified variables	Unspecified variables
Slack	[V],φ	P,Q
Generator	[V],P	Q, φ
Load(Constant power)	P,Q	[V], φ

When the power is transmitted from the generators to the loads through the transmission line the loss occurring in the transmission line is unavoidable. It is call the transmission loss and it not known in advance so there must be a provision that a generator should include the transmission loss in its generation.

This is the reason why the slack bus concept is included in the power flow analysis. The generator having largest capacity is chosen as the slack bus. Depending on the selection of the slack bus generator the power flow solution will vary.

Procedure to solve transmission system power flow

The following steps are to be followed to carry out power flow analysis. Here it is assumed that the three phase network is a balanced one.

1. Collect the specified data of generators, loads, transmission lines, and transformers. You can follow some matrix format like matpower.
2. Convert the network as per unit equivalent.
3. Formulate the equations as per kirchooff’s law.

These equations can be written in either polar form or rectangular form. In polar form no of equation are to be solved in 2(n-1)-m, wherein the rectangular form the no of equations are 2(n-1).

1. The system of non linear (polar trigonometric, rectangular quadratic) equations are to be solved by numerical methods only because so far no analytical solution is derived to sole this problem.
2.   Newton Raphson method, Decupled and gauss seidal are methods used to solve this problem. Among this N-R and Decoupled are preferred for their quadratic convergence.

Distribution Power flow.

The distribution systems are characterized by their nature of high R/X ration and radial nature (no loop).The admittance matrix is sparse and the Newton Raphson method may fail to converge for some networks. So the modeling of the problem as well as the solution methodology is different.

 The power flow equations for a radial distribution system are derived as the relationship between the specified complex bus powers and the bus voltages. If the network has n (total number) buses and ‘m’ generator buses the power flow equations are written as follows

k(i) is the set of nodes connected to node i, and P_i /Q_i denote the real/reactive power at node i. The complex non linear equations (B2) are to be solved to determine the bus voltages. So many methods ranging form sweep methods to conic programming methods available in the literature 

 to solve this problem..

.

Solving Economic Dispatch and Optimal Power Flow by GA

The economic dispatch problem is described as an act of minimizing the total fuel cost of the committed generators while satisfying the demand, network constraints. and plant.

Follow these steps.

Step II Solution Methodology

1. Choose a reference plant .For economic dispatch choose the plant with large capacity (range). In case of optimal power flow the slack bus is the reference bus.
2. For  both problems the number of  control(independent) variables to be determined in n-1.The reference plant  allocation is determined from the  constraint equations A2.

While solving the quadratic equation consider the positive solution. Check for the plant limits. If it is violating the limits allocate that particular limit.
3. Write a function file of n-1 control variables which return the total fuel cost fuel cost and the allocation. In case optimal power flow you have to use a power flow routine to determine the reference plant(slack bus) allocation.
4.The data file (fuel cost equations ,demand,loss coefficients,bus data,line data) and the gaoptions is put in a file  and the function file is run by the GA with options.


The program can be downloaded from the matlab cetral file excange

http://www.mathworks.de/matlabcentral/fileexchange/20825

Economic Load Dispatch :An Intorduction

Solving Optimization Problems Using MATLAB GA toolbox-Part 2

Now let us learn how to use the GA in command line mode.
the basic syntax to run the GA in command mode is
x = ga(@fitnessfunction,nvar,options)

fitnessfunction:Function file relating the control variables(x) with function value(F).

nvar: no of  control variables

options: Genetic algorithm parameters setting.


GA Options



PopulationType: [ 'bitstring'      | 'custom'    | {'doubleVector'} ]:-Type of the population 'custom' :  :binary,'double vector': real values

PopInitRange: [ matrix           | {[0;1]} ]: Initial value range default is [0 to 1] you can change this by giving two column vectors .For the example problem discussed for the GUI mode the limits can be like this
l=[0 0]'; u=[10 10]';

        
PopulationSize: [ positive scalar  | {20} ]-Size of the population

            
EliteCount: [ positive scalar  | {2} ]:The number best solution in each generation to be saved.

    
CrossoverFraction: [ positive scalar  | {0.8} ]: cross over value between 0 to 1

    
MigrationDirection: [ 'both'           | {'forward'} ]  The solution vector is moved in both direction.(increase or decrease)
    
MigrationInterval: [ positive scalar  | {20} ]
    
MigrationFraction: [ positive scalar  | {0.2} ]

          
Generations: [ positive scalar  | {100} ]: no of iterations
            
TimeLimit: [ positive scalar  | {Inf} ]: Time Limit for the algorithm

          
FitnessLimit: [ scalar           | {-Inf} ]:;Set value or a typical known fitness value

        
 StallGenLimit: [ positive scalar  | {50} ]: If Solution is not changing for certain number of generation stop the algorithm
        
StallTimeLimit: [ positive scalar  | {20} ]: If Solution is not changing for certain time, stop the algorithm
    
 InitialPopulation: [ matrix           | {[]} ]:You can specify initial feasible solution if already known.

  InitialScores: [ column vector    | {[]} ]: Known function values

 CreationFcn: [ function_handle  | {@gacreationuniform} ]; random function for creating initial population

 FitnessScalingFcn: [ function_handle  | @fitscalingshiftlinear  | @fitscalingprop  | :testing the fitness
                             @fitscalingtop   | {@fitscalingrank} ]
SelectionFcn: [ function_handle  | @selectionremainder    | @selectionrandom |
                            @selectionroulette | @selectiontournament   | {@selectionstochunif} ]
          
CrossoverFcn: [ function_handle  | @crossoverheuristic  | @crossoverintermediate |
@crossoversinglepoint | @crossovertwopoint | {@crossoverscattered} ]:Different crossovers
          
MutationFcn: [ function_handle  | @mutationuniform | {@mutationgaussian} ]:  Mutation function
            
HybridFcn: [ @fminsearch | @patternsearch | @fminunc | {[]} ]: After GA run another local minimizing function is run

Display: [ off | iter | diagnose | {final} ]: the results are displayed
            
OutputFcns: [ function_handle  | @gaoutputgen | {[]} ]:the output producing function
            
 PlotFcns: [ function_handle  | @gaplotbestf | @gaplotbestindiv | @gaplotdistance | @gaplotexpectation | @gaplotgeneology | @gaplotselection | @gaplotrange | @gaplotscorediversity  | @gaplotscores | @gaplotstopping  | {[]} ]: Choose the graphs you want
          
PlotInterval: [ positive scalar  | {1} ]: The plotting interval

            
Vectorized: [ 'on'  | {'off'} ]

Program in command mode.


This file  and the function file should be in the same folder which should be default.

the options can be set by the following commands

% setting the genetic algorithm parameters.
options = gaoptimset;
options = gaoptimset('PopulationSize', 50,'Generations', 500,'TimeLimit', 200,'StallTimeLimit',100,
'PlotFcns',  {@gaplotbestf,@gaplotbestindiv});
  [x ff]=ga(@ex1,2,options)

.

Solving Optimization Problems Using MATLAB GA toolbox-Part 1

The GA tool box of MATLAB is good in solving hard optimization problems. It can be run form (i) GUI (Graphical  User Interface) mode or(ii) Command line Mode.


GA A Different Introduction

Genetic Algorithm or GA is one of th basic and powerful heuristic optimization algorithms.If you read any material on this algorithm you can observe the following points.

• 1.It is working on the population dynamics.
• 2.It is searching form a set of solutions(population) .
• 3.It does not need the derivative or continuous solution search space.
• 4.It doing some probabilistic operations like Reproduction,Cross Over,Mutation on the old  population to produce a new population.

MATLAB GA Toolbox

This toolbox contains some matlab files which can do the above described actions.This can be run in two modes.The GA toolbox is written as a minimization tool. Maximization problems also can be done by converting the Maximization problem as minimization problem.

1.GUI Mode

2.Command line Mode

Steps For Solving Optimization Problem 


To use the GA toolbox you need not know anything about Ga and its dynamics. All you should know is some fundamentals about optimization or operation research and basic matlab commands.


1. Model the optimization problem as a unconstrained  minimization problem.

2. Write the unconstrained minimization problem  as matlab function file.
3. Run the matalb ga tool box in a command line mode or GUI mode.




I am giving one example in GUI .It s problem of minimizing a  two variable quadratic function subject to a linear equality constraint..Just click on the figure to enlarge the figure.

 .