Need a better guess y0 for consistent initial conditions?

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

[Torsten]

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.