Numerical workaround of false limit cycle instead of equilibrium

[strct]

Hello, I formulated my system consisting of 34 equations. It is verified in SPICE simulator, therefore it is 100% correct. Quite surprisingly it behaves numerically stable even for quite wild parameters. From this I guess jacobian is not that terrible issue here, but I don't know it surely.

First I run time domain integration to find equilibirum. Clearly it is an steady-state. But only to a given precision. In detail, there already exist some negligible oscillation. But it is of numerical nature and the system definition in fact. For me, this must be equilibrium. Because of that, the following initialization and continuation is not working.

%y0_estimation; is the last vector of my time domain integration using ODE15s
[x0,v0]=init_EP_EP(@full_DAO,y0_estimation,p,ap);


opt=contset(opt,'InitStepsize',0.0000000001);
[x_EP,v_EP,s_EP,h_EP,f_EP]=cont(@equilibrium,x0,v0,opt);

I got error that there is no convergence at x0
.
Indeed it is correct (as it appears to me from what I know). Do I understand it right?

But I was wondering If I can make some temporary workaround and set some tolerances or something to enforce EP behaviour at this particular situation.
First true, significant and actually searched for Hopf (continuing first parameter) should occur once it is more or less greater than parameter V_t value.
Background of this system is that it has just one equlibrium and the rest of the analysis is done for Hopfs and limit cycles exclusively. So continuting from equilibrium is in fact the very formal but illustrative approach. All work is then done in bifurcation of limit cycles only.
I have computed symbolically 3rd order derivatives, are they used by default then or should I allow them manually?

Full code:

init
global cds

N = 34; %number of state variables
m = 0.6;
sigma = 20;
gamma = 2;
beta = 0.08;
Vt = -1;
L1 = 40e-6;
L2 = 40e-6;
L3 = 125e-9;
L4 = 125e-9;
L5 = 125e-9;
L6 = 40e-6;
L7 = 40e-6;
L8 = 125e-9;
L9 = 125e-9;
L10 = 40e-6;
L11 = 40e-6;
L12 = (1/2 + m/2)*125e-9;
L13 = (1/2 + m/2)*125e-9;
L14 = 40e-6;
L15 = 40e-6;
L16 = 125e-9*(1 - m^2)/(2*m);
L17 = 125e-9*(1 - m^2)/(2*m);
U1 = -2;
U2 = -2;
U3 = -2;
Udd = 2.5;
Ru1 = 5;
Ru2 = 5;
Ru3 = 5;
Rudd = 2.5;
R1 = 50;
R2 = 50;
%RF loss modelling
Rc3 = 1e9;
Rc4 = 1e9;
Rc7 = 1e9;
Rc8 = 1e9;
Rc11 = 1e9;
Rc12 = 1e9;
RL3 = 1;
RL4 = 1;
RL5 = 1;
RL8 = 1;
RL9 = 1;
RL13 = 1;
RL12 = 1;
C1 = 10e-9;
C2 = 10e-9;
C3 = 50e-12;
C4 = 50e-12;
C5 = 10e-9;
C6 = 10e-9;
C7 = 50e-12;
C8 = 50e-12;
C9 = 10e-9;
C10 = 10e-9;
C11 = 50e-12;
C12 = 50e-12;
C13 = 10e-9;
C14 = m/2*50e-12;
C15 = 10e-9;
C16 = m/2*50e-12;
C17 = 10e-9;

p = [U1,U2,U3,Udd,Ru1,Ru2,Ru3,Rudd,sigma,beta,gamma,Vt,L1,L2,L3,L4,L5,L6,L7,L8,L9,L10,L11,L12,L13,L14,L15,L16,L17,C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,C11,C12,C13,C14,C15,C16,C17,R1,R2,Rc3,Rc4,Rc7,Rc8,Rc11,Rc12,RL3,RL4,RL5,RL8,RL9,RL12,RL13]';

%find equilibrium for inactive transistors
OPTIONS = [];
hls = full_DS;
[tsol,ysol] = ode15s(hls{2},[0 1000e-6],zeros(1,N),OPTIONS,U1,U2,U3,Udd,Ru1,Ru2,Ru3,Rudd,sigma,beta,gamma,Vt,L1,L2,L3,L4,L5,L6,L7,L8,L9,L10,L11,L12,L13,L14,L15,L16,L17,C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,C11,C12,C13,C14,C15,C16,C17,R1,R2,Rc3,Rc4,Rc7,Rc8,Rc11,Rc12,RL3,RL4,RL5,RL8,RL9,RL12,RL13)

figure(1)
plot(tsol,ysol)
xlabel('$t\ (s)$')
ylabel('$\mathbf{x}$')
legend
y0_estimation = ysol(end,:)';
%%
opt=contset;
opt=contset(opt,'MaxNumPoints',50);
opt=contset(opt,'MaxStepsize',1e-2);
opt=contset(opt,'Singularities',1);
opt=contset(opt,'Eigenvalues',1); %eigs of equilibirum are storen in the f array

ap=1;
opt=contset(opt,'Backward',[]);

[x0,v0]=init_EP_EP(@full_DS,y0_estimation,p,ap);

[x0(1:end-1), y0_estimation]

opt=contset(opt,'InitStepsize',0.0000000001);
[x_EP,v_EP,s_EP,h_EP,f_EP]=cont(@equilibrium,x0,v0,opt);

figure(2)
plot3(ysol(:,1),ysol(:,2),ysol(:,3))

[strct]

Here another illustration

[strct]

Also for the first parameter value -0.1 where I got periodic response, I run long integration to be sure there is a LC. Then a took like 1.5 of the period and initialize. Then I tried to continue from limic cycle but I am not getting any convergence. What can I do in that case?

%running it previsously long integration until it hit steady-state, i %just want to take 1.6 part of period
[tsol2,ysol2] = ode15s(hls{2},[0 12e-9],y0_estimation,OPTIONS,U1,U2,U3,Udd,Ru1,Ru2,Ru3,Rudd,sigma,beta,gamma,Vt,L1,L2,L3,L4,L5,L6,L7,L8,L9,L10,L11,L12,L13,L14,L15,L16,L17,C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,C11,C12,C13,C14,C15,C16,C17,R1,R2,Rc3,Rc4,Rc7,Rc8,Rc11,Rc12,RL3,RL4,RL5,RL8,RL9,RL12,RL13);
figure(3)
plot(tsol2,ysol2(:,5))
figure(4)
plot3(ysol2(:,5),ysol2(:,6),ysol2(:,8))

ap = 1;
tolerance=1e-2;
[x0,v0]=initOrbLC(@full_DS,tsol2,ysol2',p,ap,40,4,tolerance);
opt=contset;
opt=contset(opt,'MaxNumPoints',3);
%opt=contset(opt,’TSearchOrder’,0);
%opt=contset(opt,’Backward’,1);
[xlcc,vlcc,slcc,hlcc,flcc]=cont(@limitcycle,x0,v0,opt);

[Alois Steindl]

Hello,
if your simulations indicate small oscillations, then these oscillations do occur in the numerical formulation and matcont will have to take them into account.
I would bet that your system is vulnerable to round-off errors or something similar happens.
It would certainly be necessary to non-dimensionalize the model and rescale it, such that (in the ideal case) all quantities are of order 1 and also the equations contain expressions of order 1. In your parameter list there occur very small quantities.
With a poorly scaled model you will quite likely not succeed with any software package.
Good luck
Alois

[strct]

Hello Alois,

this in fact is great answear and I already started work on it. I recognised, that the problem studied previously was working perfectly for normalised values, really no problem. After renormalization, it ended up with the same error. From this I suspect that it is value dependent and I totally understand it. And yes, these values are really small, so the nondimensionalisation is the key. Or at least scaling the system.
Thank you anyway for your kind suggestion, maybe it will help someone else.