what my problem in the following matlab function.can you correct me please?
% Start of file (run_LF.m)
clear;
format short;
tole=1e-9; deg=180/pi; rad=1/deg ;
A=1.5;
% Base values
Sbase=100; Ubase=220; Zb=Ubase^2/Sbase;
% Bus data
% Number of buses
nbus=6;
%Bus 1, slack bus
U1=220/Ubase; theta1=-1.5*rad; PL1=50/Sbase; QL1=5/Sbase;
%Bus 2, PQ-bus
PG2=100/Sbase; QG2=10/Sbase; PL2=(200+0.1*A)/Sbase; QL2=20/Sbase;
%Bus 3, PQ-bus
PG3=150/Sbase; QG3=15/Sbase; PL3=300/Sbase; QL3=30/Sbase;
%Bus 4, PQ-bus
PG4=0/Sbase; QG4=0/Sbase; PL4=0/Sbase; QL4=0/Sbase; Qc=(25-0.2*A)/Sbase;
%Bus 5, PQ-bus
PG5=200/Sbase; QG5=20/Sbase; PL5=350/Sbase; QL5=(40-0.1*A)/Sbase;
%Bus 6, PQ-bus
PG6=200/Sbase; QG6=10/Sbase; PL6=100/Sbase; QL6=20/Sbase;
% Line data
r=0.0529;x=0.529;b=3.3081*10^-6;
l12=150;l13=170;l23=180;l34=100;l45=200;l56=100;
r12=r*l12;r13=r*l13;r23=r*l23;r34=r*l34;r45=r*l45;r56=r*l56;
x12=x*l12;x13=x*l13;x23=x*l23;x34=x*l34;x45=x*l45;x56=x*l56;
z12=(r12+x12*i)/Zb;b12=b*l12*Zb*j/2;
z13=(r13+x13*i)/Zb;b13=b*l13*Zb*j/2;
z23=(r23+x23*i)/Zb;b23=b*l23*Zb*j/2;
z34=(r34+x34*i)/Zb;b34=b*l34*Zb*j/2;
z45=(r45+x45*i)/Zb;b45=b*l45*Zb*j/2;
z56=(r56+x56*i)/Zb;b56=b*l56*Zb*j/2;
b21=b12;b31=b13;b32=b23;b43=b34;b54=b45;b65=b56;
% Shunt
Qshunt=Qc/Sbase;
Zshunt=Ubase^2/(Qshunt*j);
Zbshunt=Zshunt/Zb;
yshunt=i*Qshunt/1^2;
% YBUS matrix
y11=1/z12+1/z13+b12+b13; y12=-1/z12; y13=-1/z13; y14=0; y15=0; y16=0;
y21=-1/z12; y22=1/z12+1/z23+b12+b23; y23=-1/z23; y24=0; y25=0; y26=0;
y31=-1/z13; y32=-1/z23; y33=1/z13+1/z23+1/z34+b13+b23+b34; y34=-1/z34;
y35=0; y36=0;
y41=0; y42=0; y43=-1/z34; y44=1/z34+yshunt+1/z45+b34+b45; y45=-1/z45; y46=0;
y51=0; y52=0; y53=0; y54=-1/z45; y55=1/z45+1/z56+b45+b56; y56=-1/z56;
y61=0; y62=0; y63=0; y64=0; y65=-1/z56; y66=1/z56+b56;
YBUS=[ y11 y12 y13 y14 y15 y16;
y21 y22 y23 y24 y25 y26;
y31 y32 y33 y34 y35 y36;
y41 y42 y43 y44 y45 y46;
y51 y52 y53 y54 y55 y56;
y61 y62 y63 y64 y65 y66];% Define YBUS
G=real(YBUS);
B=imag(YBUS);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PGD and QGD for PU- and PQ-buses
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PGD2=PG2-PL2; % for bus 2 (PQ-bus)
PGD3=PG3-PL3; % for bus 3 (PQ-bus)
PGD4=PG4-PL4; % for bus 4 (PQ-bus)
PGD5=PG5-PL5;
PGD6=PG6-PL6;
PGD=[PGD2 ; PGD3 ; PGD4 ; PGD5 ; PGD6];
QGD2=QG2-QL2;
QGD3=QG3-QL3;
QGD4=QG4-QL4;
QGD5=QG5-QL5;
QGD6=QG6-QL6;
QGD=[QGD2 ; QGD3 ; QGD4 ; QGD5 ; QGD6];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use fsolve function in MATLAB to run load flow
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Unknown variables [theta2 theta3 theta4 U2 U3 U4]’;
% Define the initial values of the unknown variables
X0=[0 0 0 0 0 1 1 1 1 1]'; % Flat initial values
s_z=size(X0);
nx=s_z(1,1); % number of unknown variables
% The function below is used for fsolve (type "help fsolve" in MATLAB)
options_solve = optimset('Display','off','TolX',tole,'TolFun',tole);
% Parameters used for fsolve
PAR=[nx ; nbus ; U1 ; theta1];% U1 and theta1 are known (slack bus).
[X_X,FVAL,EXITFLAG,OUTPUT]=fsolve('EG2100_solve_LF',X0,options_solve,G,B,PGD,QGD,PAR);
if EXITFLAG~=1,
disp('No solution'),
EXITFLAG=EXITFLAG,
return
end,
% Solved variables X_X=[theta2 theta3 theta4 U2 U3 U4]’;
ANG=[theta1 X_X(1) X_X(2) X_X(3) X_X(4) X_X(5)];% Voltage phase angles
VOLT=[U1 X_X(6) X_X(7) X_X(8) X_X(9) X_X(10)]; % Voltage magnitudes
ANG_deg=ANG*deg;% in degrees
VOLT_kV=VOLT*Ubase; % in kV
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The generated active and reactive power at slack bus and PU-buses
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
g=-G;b=-B;
% At slack bus, bus 1
P12=g(1,2)*VOLT(1)^2-...
VOLT(1)*VOLT(2)*(g(1,2)*cos(ANG(1)-ANG(2))+b(1,2)*sin(ANG(1)-ANG(2)));
P13=g(1,3)*VOLT(1)^2-...
VOLT(1)*VOLT(3)*(g(1,3)*cos(ANG(1)-ANG(3))+b(1,3)*sin(ANG(1)-ANG(3)));
Q12=(-b12-b(1,2))*VOLT(1)^2-...
VOLT(1)*VOLT(2)*(g(1,2)*sin(ANG(1)-ANG(2))-b(1,2)*cos(ANG(1)-ANG(2)));
Q13=(-b13-b(1,3))*VOLT(1)^2-...
VOLT(1)*VOLT(3)*(g(1,3)*sin(ANG(1)-ANG(3))-b(1,3)*cos(ANG(1)-ANG(3)));
PG1=P12+P13+PL1;
QG1=Q12+Q13+QL1;
PG1_MW=PG1*Sbase;
QG1_MW=QG1*Sbase;
%%%%%%%%%
% Losses
%%%%%%%%%
P34=g(3,4)*VOLT(3)^2-...
VOLT(3)*VOLT(4)*(g(3,4)*cos(ANG(3)-ANG(4))+b(3,4)*sin(ANG(3)-ANG(4)));
PLoss_tot=(PG1+PG2+PG3+PG4+PG5+PG6)-(PL1+PL2+PL3+PL4+PL5+PL6);
PLoss_tot_MW=PLoss_tot*Sbase
PLoss_Sys1=(PG1+PG2+PG3)-(PL1+PL2+PL3+P34);% Find the power losses in
System 1 (pu)
PLoss_Sys1_MW=PLoss_Sys1*Sbase
PLoss_Sys2=(PG4+PG5+PG6+P34)-(PL4+PL5+PL6);% Find the power losses in
System 2 (MW)
PLoss_Sys2_MW=PLoss_Sys2*Sbase
% End of file (run_LF.m)
The function is in the second file named EG2100_solve_LF
% Start of file (solve_lf.m)
% This function solves g(x)=0 for x.
function [g_x]=EG2100_solve_lf(X,G,B,PGD,QGD,PAR);
nx=PAR(1); nbus=PAR(2); U1=PAR(3); theta1=PAR(4);
PGD2=PGD(1); PGD3=PGD(2); PGD4=PGD(3); PGD5=PGD(4); PGD6=PGD(5);
QGD2=QGD(1); QGD3=QGD(2); QGD4=QGD(3); QGD5=QGD(4); QGD6=QGD(5);
theta2=X(1); theta3=X(2); theta4=X(3); theta5=X(4); theta6=X(5);
U2=X(6); U3=X(7); U4=X(8); U5=X(9); U6=X(10);
ANG=[theta1 theta2 theta3 theta4 theta5 theta6]' ; VOLT=[U1 U2 U3 U4 U5
U6]';
% We have nx unknown variables, therefore the
% size of g(x) is nx by 1.
g_x=zeros(nx,1);
%Based on equation (7.43), find Pk and Qk
P2=VOLT(2)*(VOLT(1)*(G(2,1)*cos(ANG(2)-ANG(1))+B(2,1)*sin(ANG(2)-ANG(1)))...
+VOLT(3)*(G(2,3)*cos(ANG(2)-ANG(3))+B(2,3)*sin(ANG(2)-ANG(3)))...
+VOLT(2)*(G(2,2)*cos(ANG(2)-ANG(2))+B(2,2)*sin(ANG(2)-ANG(2))));
% +U2*U4*(G(2,4)*cos(theta2-theta4)+B(2,4)*sin(theta2-theta4))...
% +U2*U5*(G(2,5)*cos(theta2-theta5)+B(2,5)*sin(theta2-theta5))...
% +U2*U6*(G(2,6)*cos(theta2-theta6)+B(2,6)*sin(theta2-theta6));
P3=VOLT(3)*(VOLT(2)*(G(3,2)*cos(ANG(3)-ANG(2))+B(3,2)*sin(ANG(3)-ANG(2)))...
+VOLT(4)*(G(3,4)*cos(ANG(3)-ANG(4))+B(3,4)*sin(ANG(3)-ANG(4)))...
+VOLT(3)*(G(3,3)*cos(ANG(3)-ANG(3))+B(3,3)*sin(ANG(3)-ANG(3)))...
+VOLT(1)*(G(3,1)*cos(ANG(3)-ANG(1))+B(3,1)*sin(ANG(3)-ANG(1))));
% +U3*U6*(G(3,6)*cos(theta3-theta6)+B(3,6)*sin(theta3-theta6))...;
P4=VOLT(4)*(VOLT(3)*(G(4,3)*cos(ANG(4)-ANG(3))+B(4,3)*sin(ANG(4)-ANG(3)))...
+VOLT(5)*(G(4,5)*cos(ANG(4)-ANG(5))+B(4,5)*sin(ANG(4)-ANG(5)))...
+VOLT(4)*(G(4,4)*cos(ANG(4)-ANG(4))+B(4,4)*sin(ANG(4)-ANG(4))));
P5=VOLT(5)*(VOLT(4)*(G(5,4)*cos(ANG(5)-ANG(4))+B(5,4)*sin(ANG(5)-ANG(4)))...
+VOLT(6)*(G(5,6)*cos(ANG(5)-ANG(6))+B(5,6)*sin(ANG(5)-ANG(6)))...
+VOLT(5)*(G(5,5)*cos(ANG(5)-ANG(3))+B(5,5)*sin(ANG(5)-ANG(5))));
P6=VOLT(6)*(VOLT(5)*(G(6,5)*cos(ANG(6)-ANG(5))+B(6,5)*sin(ANG(6)-ANG(5)))...
+VOLT(6)*(G(6,6)*cos(ANG(2)-ANG(3))+B(2,3)*sin(ANG(2)-ANG(3))));
%Qk
Q2=VOLT(2)*(VOLT(1)*(G(2,1)*sin(ANG(2)-ANG(1))-B(2,1)*cos(ANG(2)-ANG(1)))...
+VOLT(3)*(G(2,3)*sin(ANG(2)-ANG(3))-B(2,3)*cos(ANG(2)-ANG(3)))...
+VOLT(2)*(G(2,2)*sin(ANG(2)-ANG(2))-B(2,2)*cos(ANG(2)-ANG(2))));
Q3=VOLT(3)*(VOLT(2)*(G(3,2)*sin(ANG(3)-ANG(2))-B(3,2)*cos(ANG(3)-ANG(2)))...
+VOLT(4)*(G(3,4)*sin(ANG(3)-ANG(4))-B(3,4)*cos(ANG(3)-ANG(4)))...
+VOLT(3)*(G(3,3)*sin(ANG(3)-ANG(3))-B(3,3)*cos(ANG(3)-ANG(3))));
Q4=VOLT(4)*(VOLT(3)*(G(4,3)*sin(ANG(4)-ANG(3))-B(4,3)*cos(ANG(4)-ANG(3)))...
+VOLT(5)*(G(4,5)*sin(ANG(4)-ANG(5))-B(4,5)*cos(ANG(4)-ANG(5)))...
+VOLT(4)*(G(4,4)*sin(ANG(4)-ANG(4))-B(4,4)*cos(ANG(4)-ANG(4))));
Q5=VOLT(5)*(VOLT(4)*(G(5,4)*sin(ANG(5)-ANG(4))-B(5,4)*cos(ANG(5)-ANG(4)))...
+VOLT(6)*(G(5,6)*sin(ANG(5)-ANG(6))-B(5,6)*cos(ANG(5)-ANG(6)))...
+VOLT(5)*(G(5,5)*sin(ANG(5)-ANG(5))-B(5,5)*cos(ANG(5)-ANG(5))));
Q6=VOLT(6)*(VOLT(5)*(G(6,5)*sin(ANG(6)-ANG(5))-B(6,5)*cos(ANG(6)-ANG(5)))...
+VOLT(6)*(G(2,3)*sin(ANG(6)-ANG(6))-B(6,6)*cos(ANG(6)-ANG(6))));
% Active power mismatch (PU- and PQ-buses)
% Bus 2
g_x(1)=P2-PGD2;
% Bus 3
g_x(2)=P3-PGD3;
% Bus 4
g_x(3)=P4-PGD4;
% Bus 5
g_x(4)=P5-PGD5;
% Bus 6
g_x(5)=P6-PGD6;
% Reactive power mismatch (PQ-buses)
% Bus 2
g_x(6)=Q2-QGD2;
% Bus 3
g_x(7)=Q3-QGD3;
% Bus 4
g_x(8)=Q4-QGD4;
% Bus 5
g_x(9)=Q5-QGD5;
% Bus 6
g_x(10)=Q6-QGD6;
% End of file (solve_lf.m)
Teshome Hambissa
Lecturer
Addis Ababa Universiy, AAiT, SECE
Addis Ababa, Ethiopia
On Tue, Apr 14, 2020 at 4:31 PM Teshome Hambissa <[email protected]>
wrote:
> Dear member,
> I need run power flow calculation for three generator and four bus
> system.can you give me the way i understand the matlab code for this problem
>
> kind regards,
>
> Teshome Hambissa
> Lecturer
> Addis Ababa Universiy, AAiT, SECE
> Addis Ababa, Ethiopia
>
>
>
> On Mon, Apr 13, 2020 at 9:58 PM Ray Daniel Zimmerman <[email protected]>
> wrote:
>
>> If I understand correctly, persistent should work if it appears at the
>> top-level within the function (as opposed to inside if-then blocks like
>> this particular case). So I think you should be able to fix this one by
>> moving that line to the first line in the function. I don’t see any other
>> persistent calls that should require changing. That is, I believe they
>> are all already at the top level of their function.
>>
>> As far as num2cell(), my first idea is to try replacing all calls with
>> a call to a new function, say mp_num2cell() that simply creates a cell
>> array of the same dimensions as the input and copies the contents element
>> by element, something like …
>>
>> function c = mp_num2cell(num)
>> c = cell(size(num));
>> for k = 1:length(num(:))
>> c{k} = num(k);
>> end
>>
>> Let me know if that works and if anything else pops up.
>>
>> Ray
>>
>>
>>
>> On Apr 13, 2020, at 9:58 AM, Lode De Herdt <[email protected]>
>> wrote:
>>
>> Hi Ray,
>>
>> Thank you for the reply. I've spent a little time trying to remove the
>> things that are unsupported from the MATPOWER code but it looks like I'm
>> opening a Pandora's box here; each error I solve just leads to the next
>> errors so it's difficult to solve without good knowledge of MATLAB and the
>> inner workings of MATPOWER.
>>
>> The first error was:
>> *PERSISTENT declarations may only appear at the top-level.*
>> *Function 'moption.m' (#57.41253.41263), line 697, column 9:*
>> *"persistent"*
>>
>> So I replaced *persistent nsc_opt* with *nsc_opt = [];* and did the same
>> for other places where 'persistent' was used as was suggested in message
>> 05192 <https://www.mail-archive.com/[email protected]/msg05192.html> of
>> this forum, but that then leads to this:
>>
>>
>> *This assignment writes a 'struct' value into a 'double' type. Code
>> generation does not support changing types through assignment. Check
>> preceding assignments or input type specifications for type mismatches.
>> Function 'mpoption.m' (#34.71098.71101), line 1531, column 5: "opt"*
>>
>> One of the other errors indicates that the *num2cell * function is not
>> supported for code generation, but who knows what else will appear once
>> that one is solved.
>>
>> Let me know if you or anyone else has any other ideas, or experience with
>> MATPOWER and code generation.
>>
>> Best regards,
>> Lode
>>
>>
>>
>> On Mon, Apr 13, 2020 at 4:36 PM Ray Daniel Zimmerman <[email protected]>
>> wrote:
>>
>>> I can’t think of another option. I’ve not personally worked with the C
>>> code generation in MATLAB, so I’m not familiar with the requirements. Do
>>> you know what specifically is used by MATPOWER that is not supported? If
>>> there are easy ways to avoid those things, it would certainly make sense to
>>> do so.
>>>
>>> Ray
>>>
>>>
>>> On Apr 13, 2020, at 7:53 AM, Lode De Herdt <[email protected]>
>>> wrote:
>>>
>>> Dear Matpower users,
>>>
>>> I am trying to run a simulation of a distribution grid on an OPAL-RT
>>> real-time digital simulator, using Matpower. The simulator requires a
>>> Simulink model, which it converts to C code in order to run the model on
>>> its hardware. Using the *MATLAB Function* block in Simulink, I can
>>> integrate a few lines of Matlab code into the model that run a power flow
>>> with Matpower.
>>>
>>> Unfortunately it seems that not all of the code used for Matpower is
>>> supported for C code generation, so the Simulink model cannot run.
>>>
>>> Other than manually re-writing the Matpower code to avoid everything
>>> that's unsupported for code generation, is there any easier way of running
>>> Matpower in Simulink? It must be compatible with code generation so using
>>> an extrinsic coder or something similar is unfortunately not an option.
>>>
>>> Thank you for your help!
>>>
>>> Best regards,
>>> Lode
>>>
>>>
>>>
>>
