Any one clear my idea to implement prohibited operation zone vector for 6 generators to get economic dispatch solution?

I am also facing the problem of implementing these Prohibited operation zone values as; Pg1-> [210 240] [350 380]; Pg2-> [90 110] [140 160]; Pg3-> [150 170] [210 240]; Pg4-> [80 90] [110 120]; Pg5-> [90 110] [140 150]; Pg6-> [75 85] [100 105]; for 6 generators using LSA. My problem is that I want to try these vector values with ramp limit coefficients. kindly help me in this to implement POZ values in LSA Code;

Answer 1

[Abhishek Kumar Singh]

Hello Asif,

To incorporate the Prohibited Operation Zones (POZs) constraints, you need to modify the Lambda Search Algorithm (LSA), which I assume you are referring to.

You can begin by defining the POZs and modify the objective functions to ensure that POZs are treated as additional constraints within this function. During the optimization process, check if the proposed generation level for each generator falls within any of its POZs.

Continue iterating, and check that all generation levels are feasible and comply with ramp rates.

Here's a highly simplified code for illustration purpose you can refer to. In a practical scenario, you'd likely use a more sophisticated optimization.

function economic_dispatch_with_POZ

% Define the demand
demand = 1000; % Total demand in MW

% Define the generator max and min outputs (in MW)
Pmax = [600, 400, 300, 200, 150, 100];
Pmin = [100, 50, 80, 50, 30, 20];

% Define the prohibited operation zones (POZs) for each generator
POZ = {
    [210, 240; 350, 380], % POZs for generator 1
    [90, 110; 140, 160],  % POZs for generator 2
    [150, 170; 210, 240], % POZs for generator 3
    [80, 90; 110, 120],   % POZs for generator 4
    [90, 110; 140, 150],  % POZs for generator 5
    [75, 85; 100, 105]    % POZs for generator 6
};

% Define ramp rate limits (in MW/h)
ramp_up_limit = [40, 30, 35, 25, 20, 15];
ramp_down_limit = [40, 30, 35, 25, 20, 15];

% Initialize lambda (cost coefficient)
lambda = 10; % Starting value for lambda
lambda_increment = 0.1; % Increment value for lambda

% Initialize generation outputs
Pg = Pmin;

% Start the lambda iteration
iteration_limit = 1000; % Limit to prevent infinite loop
iteration_count = 0;

% Start the lambda iteration
while iteration_count < iteration_limit

    iteration_count = iteration_count + 1;
    % Calculate total generation
    total_generation = sum(Pg);
    
    % Check if total generation meets demand
    if abs(total_generation - demand) < 1e-1
        % Demand is met, check for convergence of last three Pg values
        if all(abs(Pg(end-2:end) - previous_Pg(end-2:end)) < 1e-3)
            % Convergence criteria met, break the loop
            disp(['Optimal numbers found at iteration: ', num2str(iteration_count)]);
            break;
        end
    end
    
    % Adjust lambda based on whether we need more or less total generation
    if total_generation < demand
        lambda = lambda + lambda_increment;
    else
        lambda = lambda - lambda_increment;
    end
    
    % Update generation outputs based on lambda
    for i = 1:length(Pg)
        % Calculate new generation output
        new_Pg = calculate_generation_output(lambda, Pmin(i), Pmax(i));
        
        % Check ramp rate limits
        new_Pg = max(Pg(i) - ramp_down_limit(i), min(new_Pg, Pg(i) + ramp_up_limit(i)));
        
        % Check POZs and adjust generation output if necessary
        for j = 1:size(POZ{i}, 1)
            if new_Pg > POZ{i}(j, 1) && new_Pg < POZ{i}(j, 2)
                if abs(new_Pg - POZ{i}(j, 1)) < abs(new_Pg - POZ{i}(j, 2))
                    new_Pg = POZ{i}(j, 1); % Adjust to just below POZ lower bound
                else
                    new_Pg = POZ{i}(j, 2); % Adjust to just above POZ upper bound
                end
            end
        end
        
        % Update generation output
        Pg(i) = new_Pg;
    end
end

% Display the final generation outputs
disp('Final generation outputs (MW):');
disp(Pg);
end

Below is a helper function to calcaulate generation output based on lambda. Kindly note that this is just a placeholder for the actual cost function and in practical scenario, cost curves of the generators will be used. For simplicity, a linear relationship between generation output and lambda is assumed.

function Pg = calculate_generation_output(lambda, Pmin, Pmax)
    % Helper function to calculate generation output
    Pg = lambda * (Pmax - Pmin) + Pmin;
    Pg = min(max(Pg, Pmin), Pmax);
end

I hope this helps clarify the approach to implementing POZs in an economic dispatch solution.

Thanks!

[ASIF]

Thank you so much "Abhishek Kumar Singh", actually i have to use lightning search algorithm (LSA) instead of Lambda Search Algorithm(LSA).

[Abhishek Kumar Singh]

You can refer to the attached M file for reference.

[ASIF]

Thanks again "Abhishek Kumar Singh". can you help me to implement PSO for 24 hours for economic dispatch instead of hour by hour based on 3 thermal units. here is the following attached code,

model.PD=[ 140, 150, 155, 160, 165, 170, 175, 180, 210, 230, 240, 250, 240, 220, 200, 180, 170, 185, 200, 240, 225, 190, 160, 145] ;

model.plants.pmin=[37 40 50];

model.plants.pmax=[150 160 190];

model.plants.a=[0.024 0.029 0.021];

model.plants.b=[21 20.16 20.4];

model.plants.c=[1530 992 600];

model.nplant=numel(model.plants.a);

%% ModelCalculations

function out=ModelCalculations(X,model)

a=model.plants.a;

b=model.plants.b;

c=model.plants.c;

fc=a.*(X.^2)+ b.*(X)+c;

Ptotal=sum(X);

fcost=sum(fc);

end

%% Problem Definition

model=CreateModel();

CostFunction = @(X) Mycost(X,model); % Cost Function

nVar = model.nplant; % Number of Decision Variables

VarSize = [1 nVar]; % Size of Decision Variables Matrix

VarMin = [37 40 50]; % Lower Bound of Variables

VarMax = [150 160 190]; % Upper Bound of Variables with thermal only

%% PSO Parameters

MaxIt = 1000; % Maximum Number of Iterations

nPop =800; % Population Size (Swarm Size)

% PSO Parameters

wmax=0.9; % Inertia Weight

wmin=0.4;

wdamp = 0.99; % Inertia Weight Damping Ratio

c1 = 1.5; % Personal Learning Coefficient

c2 = 2.0; % Global Learning Coefficient

% If you would like to use Constriction Coefficients for PSO,

% uncomment the following block and comment the above set of parameters.

% Velocity Limits

VelMax = 0.2*(VarMax-VarMin);

VelMin = -VelMax;

%% Initialization

empty_particle.Position = [];

empty_particle.Cost = [];

empty_particle.out = [];

empty_particle.Velocity = [];

empty_particle.Best.Position = [];

empty_particle.Best.Cost = [];

empty_particle.Best.out = [];

particle = repmat(empty_particle, nPop, 1);

BestSol.Cost = inf;

for i = 1:nPop

% Initialize Position

particle(i).Position = unifrnd(VarMin, VarMax, VarSize);

% Initialize Velocity

particle(i).Velocity = zeros(VarSize);

% Evaluation

%[particle(i).Cost] = CostFunction(particle(i).Position);

[particle(i).Cost, particle(i).out] = CostFunction(particle(i).Position);

% Update Personal Best

particle(i).Best.Position = particle(i).Position;

particle(i).Best.Cost = particle(i).Cost;

particle(i).Best.out = particle(i).out;

% Update Global Best

if particle(i).Best.Cost<BestSol.Cost

BestSol = particle(i).Best;

end

end

BestCost = zeros(MaxIt, 1);

%% PSO Main Loop

%fcMatrix = zeros(MaxIt, numel(BestSol.out.fc)); % Matrix to store fc values for each iteration

for it = 1:MaxIt

w=wmin+(wmax-wmin)*(MaxIt-it)/MaxIt;

for i = 1:nPop

% Update Velocity

particle(i).Velocity = w*particle(i).Velocity ...

+c1*rand(VarSize).*(particle(i).Best.Position-particle(i).Position) ...

+c2*rand(VarSize).*(BestSol.Position-particle(i).Position);

% Apply Velocity Limits

particle(i).Velocity = max(particle(i).Velocity, VelMin);

particle(i).Velocity = min(particle(i).Velocity, VelMax);

% Update Position

particle(i).Position = particle(i).Position + particle(i).Velocity;

% Velocity Mirror Effect

IsOutside = (particle(i).Position<VarMin | particle(i).Position>VarMax);

particle(i).Velocity(IsOutside) = -particle(i).Velocity(IsOutside);

% Apply Position Limits

particle(i).Position = max(particle(i).Position, VarMin);

particle(i).Position = min(particle(i).Position, VarMax);

% Evaluation

%[particle(i).Cost] = CostFunction(particle(i).Position);

[particle(i).Cost, particle(i).out] = CostFunction(particle(i).Position);

% Update Personal Best

if particle(i).Cost<particle(i).Best.Cost

particle(i).Best.Position = particle(i).Position;

particle(i).Best.Cost = particle(i).Cost;

particle(i).Best.out = particle(i).out;

% Update Global Best

if particle(i).Best.Cost<BestSol.Cost

BestSol = particle(i).Best;

end

end

end

BestCost(it) = BestSol.Cost;

% Store fc values for each iteration

% fcMatrix(it, :) = BestSol.out.fc;

% Display the power information for each variable

disp(['Iteration ' num2str(it) ': Best Cost = ' num2str(BestCost(it))]);

w = w*wdamp;

end