Why am I getting an error 'index greater than 1' when trying to simulate a circuit based on ode15s

I am trying to simulate a circuit in a simulink model but I get an error saying that index > 1; if I change to other solvers like ode23tb or ode23s, it will have sigularity warning and the results are NaN. I have checked my code but I dont know why. Can someone help me? files are as attached.

[Meme Young]

Where is Walter Roberson, the dude who has a house at the matlab forum?

[Jan]

"an error saying that index > 1" - pleyase post a copy of the comnplete error message instead of a rough paraphrazed version. Then the readers do not have to guess, where the error occurs.

I get the error:

Error using daeic12 (line 76)
This DAE appears to be of index greater than 1.
Error in ode15s (line 310)
    [y,yp,f0,dfdy,nFE,nPD,Jfac] = daeic12(odeFcn,odeArgs,t,ICtype,Mt,y,yp0,f0,...
Error in t_ode_ieeeSSR_2nd (line 5)
[t,x] = ode15s(@ode_ieeeSSR_2nd_2VsrcNoTransTest,tspan,x0,options);

Your J is singular.

[Meme Young]

@Jan Yes I am getting exactly the same error. I am trying to simulate the lower single phase circuit in the slx file using ode(), and I have checked my code very carefully, but it still reports such an error. I really dont know why... I have searched for other answers and they cant help me.

[Walter Roberson]

https://www.mathworks.com/matlabcentral/answers/102944-what-is-the-meaning-of-this-dae-appears-to-be-of-index-greater-than-1-using-ode-solvers-for-solvin

[Meme Young]

@Walter Roberson Thank you Mr Roberson, I have also found this answer, but I still dont know how to solve this problem and fix my code

Answer 1

[Meme Young]

@Walter Roberson I have checked my matrix. The answer you posted says that the (M + lamda * dF/dy) has to be non-singular for all nonzero lamda, so I tried to construct the dF/dy and M:

M = diag([ones(1,5), zeros(1,5)]);
A = [-R12 / L12 0 0 0 0 1 / L12 0 0 0 0;
    0 -R23_1 / L23_1 0 0 0 0 1 / L23_1 0 0 0;
    0 0 -R23_2 / L23_2 0 0 0 0 1 / L23_2 0 0;
    0 0 0 0 0 0 0 0 1 / C23_2 0;
    0 0 0 0 -R34 / L34 0 0 0 0 1 / L34;
    0 0 0 1 0 0 -1 1 0 0;
    0 0 0 0 0 1 1 0 0 1;
    0 0 -1 0 0 0 0 0 1 0;
    1 -1 -1 0 0 0 0 0 0 0;
    1 0 0 0 -1 0 0 0 0 0]; % dF/dy    

But I dont know how to tell if (M + lamda * A) is non-singular for all lamda ≠ 0.

[Walter Roberson]

A singular matrix has the property that its determinant is 0, so you can take the determinant of M+lambda*A and see if there are any solutions for the root other than lambda = 0. The determinant works out as lambda^7 times a cubic polynomial in lambda. Algebra says that any polynomial in strictly real coefficients must have an even number of complex-value coefficient. A cubic in real coefficients can therefore have either one or three real roots.

The implication is that if your M and A matrices are correct, and all of those R and L and C values are real-valued, then there is at least one real-valued non-zero lambda for which the determinant is 0, and therefore the matrix is singular for at least one real-valued non-zero lambda, and possibly for three of them.

However, you did not say "real-valued", you said "all nonzero".

Suppose you let 6 out of your 9 values be 1. The determinant would then be a cubic in three variables plus lambda. Solve for lambda, and you get the three roots, which are in the three variables you did not set to 1. Now take those three lambda roots, and solve for them to be 0 simultaneously. 3 equations in 3 variables, there should be solutions. And under those particular combinations of the three constants, the cubic roots you found for the determanent in lambda, all vanish to 0 -- which would mean that there are no non-zero lambda solutions for those particular values.

Thus, if you have particular constants already determined, they probably will not satisfy the requirements that the roots of the determanent are all 0, but on the other hand you can see the theoretical possibility that you just might happening to be working with a case in which the roots did turn out all 0 even though the A matrix was non-trivial. It can happen.

[Meme Young]

@Walter Roberson Mr Roberson thank you for your answer but it is too difficult for me to understand... I am just trying to solve the equations, and I have checked carefully to make sure that they match the circuit in the slx file. Since the slx file can run, it must be a set of solvable DAE, so why will it be non-solvable using ode15s()? This is really confusing

[Walter Roberson]

baseMVA = 100; baseKV = 500;
Zbase = baseKV ^ 2 / baseMVA; 
fn = 60; wn = 2 * pi * fn;

%% line data p.u.
r12 = 0.0014; x12 = 0.03;
r23_1 = 0.0067; l23_1 = 0.0739;
r23_2 = 0.0074; l23_2 = 0.08; c23_2 = 0.55 * l23_2;
r34 = 2e-4; x34 = 0.02;

%% line data nominal
R12 = Zbase * r12; L12 = Zbase * x12 / wn;
R23_1 = r23_1 * Zbase; L23_1 = l23_1 * Zbase / wn;
R23_2 = r23_2 * Zbase; L23_2 = l23_2 * Zbase / wn; C23_2 = 1 / (c23_2 * Zbase * wn);
R34 = r34 * Zbase; L34 = x34 * Zbase / wn;

M = diag([ones(1,5), zeros(1,5)]);

A = [-R12 / L12 0 0 0 0 1 / L12 0 0 0 0;
    0 -R23_1 / L23_1 0 0 0 0 1 / L23_1 0 0 0;
    0 0 -R23_2 / L23_2 0 0 0 0 1 / L23_2 0 0;
    0 0 0 0 0 0 0 0 1 / C23_2 0;
    0 0 0 0 -R34 / L34 0 0 0 0 1 / L34;
    0 0 0 1 0 0 -1 1 0 0;
    0 0 0 0 0 1 1 0 0 1;
    0 0 -1 0 0 0 0 0 1 0;
    1 -1 -1 0 0 0 0 0 0 0;
    1 0 0 0 -1 0 0 0 0 0]; % dF/dy
syms lambda

D = det(M + lambda * A);
R = solve(D./lambda^7, lambda);

vpa(R)

ans = 

subs(M, lambda, R(3))

ans = 

That is obviously singular.

Therefore, there are three lambda values (one of them real) for which M + lambda*A is singular.

If you are correct that "the (M + lamda * dF/dy) has to be non-singular for all nonzero lamda" then you do not meet the conditions, and that would imply that you cannot reduce the DAE order.

Since the slx file can run, it must be a set of solvable DAE

Your code might possibly not match the slx model.