Need a better guess y0 for consistent initial conditions?
12 次查看(过去 30 天)
显示 更早的评论
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 个评论
回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!