clc;
clear all;
warning off;
a = [0.001562; 0.00194; 0.00482];
b = [7.92; 7.85; 7.97];
c = [561; 310; 78];
PD=850;
Pgmin = [150;100;50];
Pgmax = [600;400;200];
no_of_units=length(a);
DelP = 10;              
lambda = (PD+sum( b./(2*a)))/sum( ones(length(a), 1)./(2*a)); 
P = (lambda - b)./(2*a);
for i=1:no_of_units
    if P(i)>Pgmax(i)
        if lambda > ((2*a(i)*Pgmax(i))+b(i))
            P(i)=Pgmax(i);
        end
    elseif P(i)<Pgmin(i)
        if lambda < ((2*a(i)*Pgmax(i))+b(i))
            P(i)=Pgmax(i);
        end
    end
end
P
