Need a better guess y0 for consistent initial conditions?
12 views (last 30 days)
Show older comments
Hello,
Im trying to solve a few DAE and algebraic equations using ode15s with the help of a mass matrix.
I'm unable to understand why I'm getting this error.
What better guess for initial values should i be making? Could you please let me know if I've made a mistake?
Error using daeic12 (line 166)
Need a better guess y0 for consistent initial conditions.
Error in ode15s (line 310)
[y,yp,f0,dfdy,nFE,nPD,Jfac] = daeic12(odeFcn,odeArgs,t,ICtype,Mt,y,yp0,f0,...
Error in three_machine_with_PIcontrolLoop2 (line 153)
[t2,x2] = ode15s(@TWOAXIS_three_machine_with_PIcontrolLoop2,tspan,x0,opt);%
Calling the function to solve ODEs
Please take a look at the program below.
This is the main program:
busdata = busdata3(); line = linedata3(); bus = busdata(:,1); type = busdata(:,2);
nbus = max(bus);
pv = find(type == 2 | type == 1); npv = length(pv);
pq = find(type == 3); npq = length(pq);
Z = [ 17.361, 0, 0, 17.361, 0, 0, 0, 0, 0, 1.04, 0, 0.716, 0.1903
0, 16, 0, 0, 0, 0, 16, 0, 0, 1.028, 0.1584, 1.63, 0
0, 0, 17.065, 0, 0, 0, 0, 0, 17.065, 1.05, 0.0753, 0.85, -0
17.3611, 0, 0, 39.448, 11.6841, 10.689, 0, 0, 0, 1.0302, -0.0385, -0, 0
0, 0, 0, 11.684, 17.5252, 0, 6.0920, 0, 0, 1.0015, -0.06952, -1.25, -0.5
0, 0, 0, 10.689, 0, 16.1657, 0, 0, 5.7334, 1.0226, -0.0644, -0.9, -0.3
0, 16, 0, 0, 6.0920, 0, 35.556, 13.793, 0, 1.0331, 0.06236, 0, 0
0, 0, 0, 0, 0, 0, 13.793, 23.47, 9.852, 1.0283, 0.0103, -1.0, -0.35
0, 0, 17.065, 0, 0, 5.7334, 0, 9.852, 32.2461, 1.0510, 0.0302, -0, 0 ];
Y = Z(1:9,1:9); Vt = Z(:,10); thv = Z(:,11);
P_inj = Z(:,12); Q_inj = Z(:,13);
PL = busdata(:,7); QL = busdata(:,8);
PG = P_inj + PL; QG = Q_inj + QL;
Sg = PG+1i*QG;
Sl = PL+1i*QL;
m = length(H);
%% calculating initial values
for i = 1:m
IG(i) = (PG(i) - 1i*QG(i))./conj(Vt(i).*exp(1i*thv(i))); % net generated current at bus i
IL(i) = (-PL(i) + 1i*QL(i))./conj(Vt(i).*exp(1i*thv(i))); % net load current at bus i
I(i) = (P_inj(i) - 1i*Q_inj(i))./conj(Vt(i).*exp(1i*thv(i))); % net load current at bus i
Delta(i) = angle(Vt(i)*exp(1i*thv(i)) + (Rs(i) + 1i*Xq(i)*IG(i)));
Id(i) = real(IG(i)*exp(-1i*(Delta(i)-pi/2)));
Iq(i) = imag(IG(i)*exp(-1i*(Delta(i)-pi/2)));
Vd(i) = real(Vt(i)*exp(1i*thv(i))*exp(-1i*(Delta(i)-pi/2)));
Vq(i) = imag(Vt(i)*exp(1i*thv(i))*exp(-1i*(Delta(i)-pi/2)));
Epd(i) = (Xq(i) - Xpq(i))*Iq(i);
Epq(i) = Vq(i) + Rs(i)*Iq(i) + Xpd(i)*Id(i);
Efd(i) = Epq(i) + (Xd(i) - Xpd(i))*Id(i);
VR(i) = KE(i)*Efd(i);
Rf(i) = KF(i)*Efd(i)/TF(i);
Vref(i) = Vt(i) + VR(i)/KA(i);
Tm(i) = Epd(i)*Id(i) + Epq(i)*Iq(i) + (Xpq(i) - Xpd(i))*Id(i)*Iq(i);
Psv(i) = Tm(i);
PC(i) = Psv(i);
w(i) = ws;
F_i(i) = 0;
F_coi(i) = 0;
end
%% constructing the mass matrix that is needed for solving the DAE
MM_machine = [1 1 1 1 1 1 1 1 1 1 1]'; % Mass matrix addition for 1 machine
MM_machine_alegraic = [0 0]'; %for algrebraic equations
MM_bus = [0 0]'; % Mass matrix addition for 1 bus -- one zero each of the...
% ... matrix for the two power balance equations
MM = MM_machine; % Generates Mass Matrix for m-machine n-bus system
m=length(H);
if m > 1 for i = 2:m MM = vertcat(MM,MM_machine); end end
m=length(H);
if m > 1 for i = 1:m MM = vertcat(MM,MM_machine_alegraic); end end
if nbus > 1 for i = 1:nbus MM = vertcat(MM,MM_bus); end end
%% Calling the function file to solve the system of DAE
st = 0.0164; %step
tf = 60; %final time
tspan = [ 0 tf ]; %time span
M = diag( MM ); % Mass matrix (forms a diagonal matrix)
opt = odeset( 'Mass',M ,'MaxStep',st); % ODE solver option
x0 = [Epq Epd Delta w Efd Rf VR Tm Psv F_i F_coi Id Iq Vt' thv' ]'; % Initial Values or S.S Values
[t2,x2] = ode15s(@TWOAXIS_three_machine_with_PIcontrolLoop2,tspan,x0,opt);% Calling the function to solve ODEs
and this is the DAE function file:
function xd = TWOAXIS_three_machine_with_PIcontrolLoop2(t,x)
global tf Xd Xpd Xq Xpq ws Rs Tpd0 Tpq0 H TFW KE TE KF TF KA Vref TCH TA PC...
RD Tsv Y Sl nbus type npv npq KP KI
xd = zeros(15*npv+2*npq,1);
%% State space variables:
Epq = x(1:npv); Epd = x(npv+1:2*npv); Del = x(2*npv+1:3*npv); w = x(3*npv+1:4*npv);
Efd = x(4*npv+1:5*npv); Rf = x(5*npv+1:6*npv); VR = x(6*npv+1:7*npv); Tm = x(7*npv+1:8*npv);
PSV = x(8*npv+1:9*npv); F_I = x(9*npv+1:10*npv); F_coi = x(10*npv+1:11*npv); Id = x(11*npv+1:12*npv);
Iq = x(12*npv+1:13*npv); Vt = x(13*npv+1:14*npv+npq); thv = x(14*npv+npq+1:14*npv+npq+nbus);
Vref2 = 1;
PC2 = 1;
%% Center Of Inertia
for ii=1:npv
w1(ii)=H(ii)*(w(ii)-ws);
end
w_coi = sum(w1)/sum(H);
%% DAE to be solved
for i=1:npv
if type(i)==1 || type(i)==2 %%%%%%
xd(i,1) = (-Epq(i)-((Xd(i)-Xpd(i))*Id(i))+Efd(i))/Tpd0(i); %
xd(i,2) = (-Epd(i)+((Xq(i)-Xpq(i))*Iq(i)))/Tpq0(i); %
xd(i,3) = (w(i)-ws); %
xd(i,4) = (Tm(i)-(Epd(i)*Id(i))-(Epq(i)*Iq(i))-((Xpq(i)-Xpd(i))*Id(i)*Iq(i))...
-TFW(i)-0.1*(w(i)-ws))*(ws/(2*H(i))); %
xd(i,5) = (-(KE(i))*Efd(i)+VR(i))/TE(i); % D
xd(i,6) = (-Rf(i)+KF(i)*Efd(i)/TF(i))/TF(i); % A
xd(i,7) = (-VR(i)+KA(i)*Rf(i)-KA(i)*KF(i)*Efd(i)/TF(i)+KA(i)*(Vref2(i)-Vt(i)))/TA(i); % E
xd(i,8) = (-Tm(i)+PSV(i))/TCH(i); %
xd(i,9) = (-PSV(i) + PC2(i) - (1/RD(i))*(w(i)/ws-1) - F_coi(i))/Tsv(i); %
xd(i,10) = KI(i)*w_coi; %
xd(i,11) = (KP(i)*w_coi) + F_I(i); %%%%%%
xd(i,12) = Epd(i)-Vt(i)*sin(Del(i)-thv(i))-Rs(i)*Id(i)+Xpq(i)*Iq(i); % Algebraic
xd(i,13) = Epq(i)-Vt(i)*cos(Del(i)-thv(i))-Rs(i)*Iq(i)-Xpd(i)*Id(i); % Equations
end
end
%% power balance equations (also Algebraic Equations)
for i =1:nbus
P(i,1) = 0;
Q(i,1) = 0;
if type(i)==1
for k=1:nbus
P(i,1) = P(i,1)+Vt(i)*Vt(k)*abs(Y(i,k))*cos(thv(i)-thv(k)-angle(Y(i,k)));
Q(i,1) = Q(i,1)+Vt(i)*Vt(k)*abs(Y(i,k))*sin(thv(i)-thv(k)-angle(Y(i,k)));
end
xd(i,14) = Id(i)*Vt(i)*sin(Del(i)-thv(i))+Iq(i)*Vt(i)*cos(Del(i)-thv(i))-real(Sl(i))-P(i);
xd(i,15) = Id(i)*Vt(i)*cos(Del(i)-thv(i))-Iq(i)*Vt(i)*sin(Del(i)-thv(i))-imag(Sl(i))-Q(i);
end
if type(i)==2
for k=1:nbus
P(i,1) = P(i,1)+Vt(i)*Vt(k)*abs(Y(i,k))*cos(thv(i)-thv(k)-angle(Y(i,k)));
Q(i,1) = Q(i,1)+Vt(i)*Vt(k)*abs(Y(i,k))*sin(thv(i)-thv(k)-angle(Y(i,k)));
end
xd(i,14) = Id(i)*Vt(i)*sin(Del(i)-thv(i))+Iq(i)*Vt(i)*cos(Del(i)-thv(i))-real(Sl(i))-P(i);
xd(i,15) = Id(i)*Vt(i)*cos(Del(i)-thv(i))-Iq(i)*Vt(i)*sin(Del(i)-thv(i))-imag(Sl(i))-Q(i);
end
if type(i)==3
for k=1:nbus
P(i,1) = P(i,1)+Vt(i)*Vt(k)*abs(Y(i,k))*cos(thv(i)-thv(k)-angle(Y(i,k)));
Q(i,1) = Q(i,1)+Vt(i)*Vt(k)*abs(Y(i,k))*sin(thv(i)-thv(k)-angle(Y(i,k)));
end
xd(i,14) = real(Sl(i))+P(i);
xd(i,15) = imag(Sl(i))+Q(i);
end
end
1 Comment
Answers (0)
See Also
Categories
Find more on Linear Algebra in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!