Need a better guess y0 for consistent initial conditions?

11 次查看（过去 30 天）
Joseph Nikhil Raj Koppula
编辑： 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 Center 和 File 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