code of load flow equation four buses

clc
clear all
cg =[0.86 0 0 0.6];
Pg = [0 0 0 318];
Pd = [50 170 200 80];
Qd = [30.9 105.35 123.94 49.58];
Y = [9.934-1i*44.925, -3.815+1i*19.078, -5.169+1i*25.847, 0
-3.815+1i*19.078, 8.984 – 1i*44.925, 0, -5.169+1i*25.847
-5.169+1i*25.847, 0, 8.192 – 1i*40.965, -3.023 + 1i*15.118
0, -5.169+1i*25.847, -3.023 + 1i*15.118, 8.192 – 1i*40.965];
Pg=sdpvar(1,4);
Qg = sdpvar(1,4);
objective= sum(cg.*Qg)
Qgmin=[-100 0 0 -100];
Qgmax =[100 0 0 100];
Pgmin =[0 0 0 0];
Pgmax =[318 0 0 318];
vang = sdpvar(1,4);
vmag = sdpvar(1,4);
vmin =0.9;
vmax =1.1
Constraints = [];

[Costraints,Pg(1,1)-Pd(1,1)==vmag(1,1)*[vmag(1,1)*Ymag(1,1)*cos(Yang(1,1))+…
vmag(1,2)*Ymag(2,1)*cos(Yang(2,1) + vang(1,1) – vang(1,2))+…
vmag(1,3)*Ymag(3,1)*cos(Yang(3,1) + vang(1,1) – vang(1,3))+…
vmag(1,4)*Ymag(4,1)*cos(Yang(4,1) + vang(1,1) – vang(4))];
[Costraints,Qg(1,1)-Qd(1,1)==(-vmag)(1,1)*[vmag(1,1)*Ymag(1,1)*sin(Yang(1,1))+…
vmag(1,2)*Ymag(2,1)*sin(Yang(2,1) + vang(1,1) – vang(1,2))+…
vmag(1,3)*Ymag(3,1)*sin(Yang(3,1) + vang(1,1) – vang(1,3))+…
vmag(1,4)*Ymag(4,1)*sin(Yang(4,1) + vang(1,1) – vang(4))];

% Constraints = [Constraints,sum(Pg)==sum(Pd)];
% Constraints=[Constraints,Pg(i)-Pd(i)=P(i)]
% Constraints=[Constraints,Qg(i)-Qd(i)=Q(i)]
Constraints = [Constraints,0<=Pg<=318];
Constraints = [Constraints,-100<=Qg<=100];
Constraints = [Constraints, vmin <= vmag];
Constraints = [Constraints, vmag <= vmax];
options = sdpsettings(‘solver’,’gurobi’);
sol = optimize(Constraints,objective,options);
value(objective)
reactive_power=value(Qg)
active_power=value(Pg)
voltage_mag=value(vmag)
voltage_ang=value(vang)clc
clear all
cg =[0.86 0 0 0.6];
Pg = [0 0 0 318];
Pd = [50 170 200 80];
Qd = [30.9 105.35 123.94 49.58];
Y = [9.934-1i*44.925, -3.815+1i*19.078, -5.169+1i*25.847, 0
-3.815+1i*19.078, 8.984 – 1i*44.925, 0, -5.169+1i*25.847
-5.169+1i*25.847, 0, 8.192 – 1i*40.965, -3.023 + 1i*15.118
0, -5.169+1i*25.847, -3.023 + 1i*15.118, 8.192 – 1i*40.965];
Pg=sdpvar(1,4);
Qg = sdpvar(1,4);
objective= sum(cg.*Qg)
Qgmin=[-100 0 0 -100];
Qgmax =[100 0 0 100];
Pgmin =[0 0 0 0];
Pgmax =[318 0 0 318];
vang = sdpvar(1,4);
vmag = sdpvar(1,4);
vmin =0.9;
vmax =1.1
Constraints = [];

% Constraints = [Constraints,sum(Pg)==sum(Pd)];
% Constraints=[Constraints,Pg(i)-Pd(i)=P(i)]
% Constraints=[Constraints,Qg(i)-Qd(i)=Q(i)]
Constraints = [Constraints,0<=Pg<=318];
Constraints = [Constraints,-100<=Qg<=100];
Constraints = [Constraints, vmin <= vmag];
Constraints = [Constraints, vmag <= vmax];
options = sdpsettings(‘solver’,’gurobi’);
sol = optimize(Constraints,objective,options);
value(objective)
reactive_power=value(Qg)
active_power=value(Pg)
voltage_mag=value(vmag)
voltage_ang=value(vang) clc
clear all
cg =[0.86 0 0 0.6];
Pg = [0 0 0 318];
Pd = [50 170 200 80];
Qd = [30.9 105.35 123.94 49.58];
Y = [9.934-1i*44.925, -3.815+1i*19.078, -5.169+1i*25.847, 0
-3.815+1i*19.078, 8.984 – 1i*44.925, 0, -5.169+1i*25.847
-5.169+1i*25.847, 0, 8.192 – 1i*40.965, -3.023 + 1i*15.118
0, -5.169+1i*25.847, -3.023 + 1i*15.118, 8.192 – 1i*40.965];
Pg=sdpvar(1,4);
Qg = sdpvar(1,4);
objective= sum(cg.*Qg)
Qgmin=[-100 0 0 -100];
Qgmax =[100 0 0 100];
Pgmin =[0 0 0 0];
Pgmax =[318 0 0 318];
vang = sdpvar(1,4);
vmag = sdpvar(1,4);
vmin =0.9;
vmax =1.1
Constraints = [];

Cart

Cart