Working with precision with integrals

Good morning!

I am working on a dynamic control problem with state space matrices.

I have generated a script where I describe these matrices. Although apparently they are correct, the result of some integrals are not correct.

I am afraid that it is a precision problem since these matrices, in themselves, are singular and have no inverse.

Specifically, the invalid results are those of Vc1(line 76), Vc2(line 77) and Vc3(line 78).

clear all
close all
clc

format longEng
%%
%Input the data of the converter
VH = 200;
VL = 70;

Llk = 2*10^-6;
n = 2/5;
fsw = 500*10^3;
Ts = 1/(2*fsw);  %Half-period, because DutiFi was calculated for half period.

Pout = 1080;
Ilow = Pout/VL;
Ro =  VL^2/Pout;
Rl = 100e-3;

%Required Duty is calculated
DutyFi = -(sqrt( VH*Ts*(VH*Ts - 16*Llk*Ilow*n) )-VH*Ts)/(4*VH*Ts);
phi = DutyFi*Ts;
Ts=2*Ts; %Undo the fix I had made earlier

C1=47*10^-6; 
C2=C1;
C3=47*10^-6;
C4=C3;

D1=DutyFi;
D2=0.5-DutyFi;
t1=D2*Ts;
t2=D1*Ts;
t3=D2*Ts;
t4=D1*Ts;

%%
%State-Space Matrices

%%%State-Space Matrices defined as:
%%% X = [Il, Vc1, Vc2, Vc3]
%%% V = [VH]

V = VH;
Cs = eye(4);

A1 = [-Rl/Llk 0 -1/(n*Llk) -1/Llk; 0 -1/(Ro*C1) -1/(Ro*C1) 0; 1/(n*C2) -1/(Ro*C2) -1/(Ro*C2) 0 ; 1/(C4+C3) 0 0 0];
B1 = [1/Llk; 0; 0; 0];
A2 = [-Rl/Llk 0 -1/(n*Llk) -1/Llk; 0 -1/(Ro*C1) -1/(Ro*C1) 0; 1/(n*C2) -1/(Ro*C2) -1/(Ro*C2) 0; 1/(C3+C4) 0 0 0];
B2 = [0; 0; 0; 0];
A3 = [-Rl/Llk 1/(n*Llk) 0 -1/Llk; -1/(n*C1) -1/(Ro*C1) -1/(Ro*C1) 0; 0 -1/(Ro*C2) -1/(Ro*C2) 0; 1/(C3+C4) 0 0 0];
B3 = [0; 0; 0; 0];
A4 = [-Rl/Llk 1/(n*Llk) 0 -1/Llk; -1/(n*C1) -1/(Ro*C1) -1/(Ro*C1) 0;  0 -1/(Ro*C2) -1/(Ro*C2) 0; 1/(C3+C4) 0 0 0];
B4 = [1/Llk; 0; 0; 0];

Xt4=[0; 0; 0; 0]; %%Start from repose

psi1= integral(@(t)expm(A1*t),0,t1, 'ArrayValued',1);
psi2= integral(@(t)expm(A2*t),0,t2, 'ArrayValued',1);
psi3= integral(@(t)expm(A3*t),0,t3, 'ArrayValued',1);
psi4= integral(@(t)expm(A4*t),0,t4, 'ArrayValued',1);

for g=1:20000 %If we go through all the states and repeat enough times we will reach the stationary state.
    Xt1=expm(A1*t1)*Xt4+psi1*B1*V;
    Xt2=expm(A2*t2)*Xt1+psi2*B2*V;
    Xt3=expm(A3*t3)*Xt2+psi3*B3*V;
    Xt4=expm(A4*t4)*Xt3+psi4*B4*V;
end

%Plot the waveform of each state variable to check that it is working properly.
IL=[ Xt4(1), Xt1(1), Xt2(1), Xt3(1), Xt4(1)];
Vout=[ Xt4(3)+Xt4(2), Xt1(3)+Xt1(2), Xt2(3)+Xt2(2), Xt3(3)+Xt3(2), Xt4(3)+Xt4(2)];
Vc1=[ Xt4(2), Xt1(2), Xt2(2), Xt3(2), Xt4(2)]; %Should give a result of ~35
Vc2=[ Xt4(3), Xt1(3), Xt2(3), Xt3(3), Xt4(3)];  %Should give a result of ~35
Vc3=[ 2*Xt4(4), 2*Xt1(4), 2*Xt2(4), 2*Xt3(4), 2*Xt4(4)];    %%Should give a result of ~100

t0=0;
time=[t0 t0+t1 t0+t1+t2 t0+t1+t2+t3 t0+t1+t2+t3+t4];

figure(1)
subplot(411)
plot(time, IL,'o-', 'LineWidth',2)
title('Inductor current')
xlabel('Time (Seconds)')
ylabel('Ampers(A)')
grid on
grid minor
subplot(412)
plot(time, Vc1,'o-', 'LineWidth',2)
title('Vc1')
xlabel('Time (Seconds)')
ylabel('Volts(V)')
grid on
grid minor
subplot(413)
plot(time, Vc2,'o-', 'LineWidth',2)
title('Vc2')
xlabel('Time (Seconds)')
ylabel('Volts(V)')
grid on
grid minor
subplot(414)
plot(time, Vc3,'o-', 'LineWidth',2)
title('Vc3')
xlabel('Time (Seconds)')
ylabel('Volts(V)')
grid on
grid minor

Any ideas?

Thank you very much!

Answer 1

[Walter Roberson]

I adapted your code to use symbolic computation. I used the default 32 digits for routine computation, and I used relative and absolute tolerance 1e-20 for the integrations -- so 4 more or digits more precise than what you would expect from integral()

The results for Vc1, Vc2, Vc3, were the same as your numeric results, to 12 decimal places.

You do not appear to have a precision problem.

Answer 2

[SERGIO TOLEDANO DUS]

Thank you very much for your help.

The reason for posting was that the matrices we integrate are not invertible. Therefore I get the message that it is not possible to work with precission.

I think it is related, because the same methodology with other exercises gives good results...

[Walter Roberson]

What operation are you doing that gives you the message about not invertible?

If you are talking about code that you did not post in which you hypothetically tried to find the steady state analytically, but were stopped because A1, A2, A3 A4 are not invertible, then the reason those are not invertible is because of the way you constructed them.

syms Rl Llk n Ro C1 C2 C3 C4
A1 = [-Rl/Llk 0 -1/(n*Llk) -1/Llk; 0 -1/(Ro*C1) -1/(Ro*C1) 0; 1/(n*C2) -1/(Ro*C2) -1/(Ro*C2) 0 ; 1/(C4+C3) 0 0 0]

A1 = 

rank(A1)

ans = 3

Notice that the second and third entry of the second row are the same as the second and third entry of the third row, so subtracting row two from row three gives a result which is zero except in the first column. But the fourth row is also zero except in the first column. Your fourth row is a linear combination of the second row and third row, so A1 is singular no matter how precisely you calculate.