From: Gwen Farrow on
Hi everyone, I'm new here, and a begginer at matlab; I am working in a project using the ZAD technique to control a boost converter. I've already figured out what are the equations for my problem, and I have put them into a function that calculates what i need; this function is presented below:

function X=switchinboostzad(xi,A1,A2,B,K,T,gamma,x1ref,x2ref)
Xref=[x1ref; x2ref],
s=K*(xi-Xref);
s1=K*(A1*xi+B);
s2=K*(A2*xi+B);
d=((2*s+T*s2)/(s2-s1))/T;
if 0<d & d<1
phi1=expm(A1*d/2*T);
psi1=B*d/2*T;
Y0=phi1*xi+psi1;
phi2=expm(A2*(T-d*T));
psi2=A2\(expm(A2*(T-d*T))-eye(2))*B;
Z0=phi2*Y0+psi2;
X=phi1*Z0+psi1;
elseif d<=0
d=0;
phi2=expm(A2*T);
psi2=A2\(expm(A2*T)-eye(2))*B;
X=phi2*xi+psi2;
elseif d>=1
d=1;
phi1=expm(A1*T);
psi1=B*T;
X=phi1*xi+psi1;
end

So, starting from an xi point (as a column vector), This function gives the value of X after passing through the control action: computes a duty cycle d for the pwm centered signal, and depending of this value computes the next x after a period of the pwmc signal has passed, using what is called a Poincare Map.
So, this is f(xi)=X
The other parameters such as A1,gamma etc., are constants:

gamma=.35; T=0.18; k2=-1; x1ref=2.5; x2ref=(x1ref)^(2)*gamma;
A1= [-gamma 0; 0 0]; B=[0;1]; A2= [-gamma 1;-1 0];
k1 can take any value between 0.1 and 1.

The value I took for xi was [x1ref;x2ref], and X is delivered in a column vector form too.

Well, I needed to find out what was the value for xi, that was equal to X, i.e. I need to find the roots for f(xi)=xi i.e. X=xi ; X-xi=0 or xi-X=0 and to do that, I did the next change to the function above:

function out=switchinboostzad(xi,A1,A2,B,K,T,gamma,x1ref,x2ref)
Xref=[x1ref; x2ref],
s=K*(xi-Xref);
s1=K*(A1*xi+B);
s2=K*(A2*xi+B);
d=((2*s+T*s2)/(s2-s1))/T;
if 0<d & d<1
phi1=expm(A1*d/2*T);
psi1=B*d/2*T;
Y0=phi1*xi+psi1;
phi2=expm(A2*(T-d*T));
psi2=A2\(expm(A2*(T-d*T))-eye(2))*B;
Z0=phi2*Y0+psi2;
X=phi1*Z0+psi1;
elseif d<=0
d=0;
phi2=expm(A2*T);
psi2=A2\(expm(A2*T)-eye(2))*B;
X=phi2*xi+psi2;
elseif d>=1
d=1;
phi1=expm(A1*T);
psi1=B*T;
X=phi1*xi+psi1;
end
out=[xi-X];
It is a system of nonlinear equations with 2 variables, so I used the fsolve function to find the zeros of this equation (xi-X=0, where X is a function of xi) from an starting point xi:
Xref=[x1ref; x2ref];
xi=Xref;
[x0,fval,exitflag,output,jacobian]=fsolve(@switchinboostzad,xi,[],A1,A2,B,K,T,gamma,x1ref,x2ref);
I get the next result for x0:
2.4993
2.1873
It seems to work.
Now, I need to know the solution for the second iteration of this function:
f(f(xi))=xi
And the third iteration:
f(f(f(xi)))=xi
and the fourth:
f(f(f(f(xi))))=xi
and the nth:
f(f(f(f(.....f(xi)....))))=xi n-times

(note: this iterations are not the same as the iterations of fsolve.)

And I just dont know how to do it.
I tried to do the next change to the function above:

function [out]=switchinboostzad(xi,A1,A2,B,K,T,gamma,x1ref,x2ref,numo)
Xref=[x1ref; x2ref];
for i=1:numo
if i==1
s=K*(xi-Xref);
s1=K*(A1*xi+B);
s2=K*(A2*xi+B);
d=((2*s+T*s2)/(s2-s1))/T;
if 0<d & d<1
phi1=expm(A1*d/2*T);
psi1=B*d/2*T;
Y0=phi1*xi+psi1;
phi2=expm(A2*(T-d*T));
psi2=A2\(expm(A2*(T-d*T))-eye(2))*B;
Z0=phi2*Y0+psi2;
X=phi1*Z0+psi1;
elseif d<=0
d=0;
phi2=expm(A2*T);
psi2=A2\(expm(A2*T)-eye(2))*B;
X=phi2*xi+psi2;
elseif d>=1
d=1;
phi1=expm(A1*T);
psi1=B*T;
X=phi1*xi+psi1;
end
else
s=K*(X-Xref);
s1=K*(A1*X+B);
s2=K*(A2*X+B);
d=((2*s+T*s2)/(s2-s1))/T;
if 0<d & d<1
phi1=expm(A1*d/2*T);
psi1=B*d/2*T;
Y0=phi1*X+psi1;
phi2=expm(A2*(T-d*T));
psi2=A2\(expm(A2*(T-d*T))-eye(2))*B;
Z0=phi2*Y0+psi2;
X=phi1*Z0+psi1;
elseif d<=0
d=0;
phi2=expm(A2*T);
psi2=A2\(expm(A2*T)-eye(2))*B;
X=phi2*X+psi2;
elseif d>=1
d=1;
phi1=expm(A1*T);
psi1=B*T;
X=phi1*X+psi1;
end
end
end
out=[xi-X];

where numo is the number of iterations of the function, and in the script, it goes like this:
options = optimset('Jacobian','off');
[x0,F,exitflag,output,JAC]=fsolve(@switchinboostzad,xi,options,A1,A2,B,K,T,gamma,x1ref,x2ref,numo);

but it gives the same outcome for x0 for any value of numo, although it is a bit different with format long:

for numo=1 x0=
2.499308123908802
2.187384146706227
for numo=2 x0=
2.499308014143268
2.187384088662747
for numo=3 x0=
2.499308036596760
2.187384135574587
for numo=4 x0=
2.499308026058130
2.187384092188981
for numo=10 x0=
2.499308049304966
2.187384099869770
And I'm getting very simmilar results.
The questions are: ¿Am I using this function fsolve the wrong way?
¿It is ok to iterate the function switchinboostzad this way (with a for)?
¿What happens if I want to compute the Jacobian for the Poincare map of this function and his iterations; it is ok to do it with the default function in matlab?
If I should not use fsolve this way, or cannot iterate the poincare map the way I have done it, ¿Is there any way to use fsolve to achieve this? or ¿What function or process or algorithm should I use?
¿if I use 'TolFun',1e-17 in the options of fsolve, is there a chance to get more suitable results?
I appreciate any help that anyone can give to me.