필터 지우기
필터 지우기

Need a better guess y0 for consistent initial conditions?

조회 수: 14 (최근 30일)
Joseph Nikhil Raj Koppula
Joseph Nikhil Raj Koppula 2019년 4월 2일
편집: Torsten 2019년 4월 2일
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
Torsten
Torsten 2019년 4월 2일
편집: Torsten 2019년 4월 2일
Insert the initial values for the differential variables into the algebraic equations and solve the algebraic equations for the algebraic variables. This will give you a vector of initial values which is consistent with your DAE system.

댓글을 달려면 로그인하십시오.

답변 (0개)

카테고리

Help CenterFile Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기

제품


릴리스

R2018b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by