calculate the Γ matrix in MATLAB from Φ Matrix - State space equation

Could you please advise me if there is a function to calculate the Γ matrix in MATLAB from Φ Matrix?

[Torsten]

If you don't explain what the matrices are, it's a question for insiders and unlikely that you get an answer.

[Rajesh Kanesan]

Thanks for your response. I am a beginner in the use of the MATLAB application. I have completed my assignment

Discretise the system into a discrete-time system in the form below.

X (k + 1) = Phi X (k) + Gammau(k)

y = CX (k)

I have manually calculated as per the steps below but would appreciate it if I could validate my answers using MATLAB.

Step 1 - Find Transition Matrix.

Step 2 - Find Gamma

A = [0 1; -1 -1]

b = [0;1]

I = eye(2);

syms s;

TransitionMatrix = (s*I-A)^-1

phi = ilaplace(TransitionMatrix)

A = [0 1; -1 -1]

I = eye(2);

syms s;

TransitionMatrix = (s*I-A)^-1

phi = ilaplace(TransitionMatrix)

[Sam Chak]

@Rajesh Kanesan

I don't know that you can do a discrete-time system with the complex variable s.

Please show your working in pencil and paper so that it can be validated in MATLAB.

[Paul]

So the goal is to compute phi and gamma from A and B?

Why is phi computed twice in the same way in the code above, but the code doesn't compute gamma.

Is a symbolic solution desired? Or a numerical solution for a specified value of the sampling period, T?

[Rajesh Kanesan]

Thanks for your response . Symbolic should be sufficient with the T .

Answer 1

[Paul]

Let's see what we have so far:

A = [0 1; -1 -1];
b = [0;1];
I = eye(2);
syms s;
LaplaceTransitionMatrix = (s*I-A)^-1

LaplaceTransitionMatrix = 

phi = ilaplace(LaplaceTransitionMatrix)

phi = 

Note that phi is a function of t, because that's the default variable for ilaplace(). But we can replace that with T

syms t T
phi = subs(phi,t,T)

phi = 

Now that we have phi, do you have a formula for gamma? If so, the code you're trying to use to implement that formula.

[Rajesh Kanesan]

Thanks for your help Paul. The below screenshot shows how i calculated manually for a similar question. all the T are replaced with q and calculated as below.

[Paul]

The equation for gamma is correct, I believe, but its implementation seems incorrect. Repeating from above to compute phi

A = [0 1; -1 -1];
b = [0;1];
I = eye(2);
syms s
LaplaceTransitionMatrix = (s*I-A)^-1;
phi = ilaplace(LaplaceTransitionMatrix);
syms t T
phi = subs(phi,t,T);

Express phi to three decimal places.

vpa(phi,3)

ans = 

which should be the matrix in the integrand for gamma (with T replaced with q). So I'm not sure where the terms in exp(-100*q) and exp(-1.71*q) come from. Maybe that was for a different A matrix?

[Rajesh Kanesan]

Thanks Paul , that was for the different Matirx . I simulated the system as below in MATLAB for that question . I wasn't able to validate the Gamma matrix in the MatLab to verify my calculation.

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% %%%%%%%%%%%%%%% Assignment 3 - Question 3 Matlab Code %%%%%%%%%%%%%%%%%%

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% File Name : Assignment1.m

%% Authors Name : Rajesh Kanesan, from MEPR, USQ

%% Email Address : [e-mail]

%% Environment : Matlab R2015b

%% Log of Change : 2022 S1, Initial Version

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Simulation of the discrete system in state space equations

% State Vector: X(k+1) = Phi*X(k) + Gamma*u(k);

% | | |

% current previous previous

% system output: y(k) = C*X(k);

clc;

clear all;

%Constants

T = 0.1;

a = (-0.017*exp(-100*T)+1.017*exp(-1.71*T));

b = (-0.010*exp(-100*T)+0.010*exp(-1.71*T));

c = (1.734*exp(-100*T)-1.734*exp(-1.71*T));

d = (1.017*exp(-100*T)-0.017*exp(-1.71*T));

e = [0.0001*exp(-100*T)-0.0058*exp(-1.71*T)+0.0057];

f = [-0.0102*exp(-100*T)+0.0099*exp(-1.71*T)+0.0002];

Phi = [a b;c d]; %Plant Matrix

Gamma = [e;f]; %Driving Matrix

C = [20.83 0]; %Connection Matrix

UnitInputs = ones(1,100); %100 Unit Step Inputs

Outputs = zeros(1,100); %100 Corresponding Outputs

StateVectorX = cell(1,100); %100 Corresponding State Vector

for index = 1:100

if index == 1 % Inital start without previous information

StateVectorX{1,index} = [0;0];

Outputs(1,index) = C*StateVectorX{1,index};

else

StateVectorX{1,index} = Phi*StateVectorX{1,index-1}+Gamma*UnitInputs(1,index-1);

Outputs(1,index) = C*StateVectorX{1,index};

end

end

plot([0:99],Outputs,'.-b');% plot output

grid on;

title('Simulation of the Discrete System - Question 3'); % plot title

legend(['T Value = ',num2str(T)]); % plot legend

xlabel('t steps'), ylabel('Output y(k)'); % plot axis labels

[Paul]

Ok. Let's see how we can solve the problem a couple of ways in Matlab

A = [0 1; -1 -1];
B = [0;1];
I = eye(2);
syms s
LaplaceTransitionMatrix = (s*I-A)^-1;
phi = ilaplace(LaplaceTransitionMatrix);
syms t T
phi = subs(phi,t,T);

Compute gamma via its defining integral

syms tau
gamma1 = simplify(int(subs(phi,T,tau)*B,tau,0,T))

gamma1 = 

With A invertible, compute gamma via the very cool equation provide by @Sam Chak

gamma2 = simplify(A\(phi - I)*B)

gamma2 = 

We can also compute phi and gamma simultaneously

temp = expm([A*T B*T;zeros(1,3)]);
phi3 = simplify(rewrite(temp(1:2,1:2),'sincos'),100)

phi3 = 

gamma3 = simplify(expand(rewrite(temp(1:2,3),'sincos')),100,'Criterion','PreferReal')

gamma3 = 

I don't know why it's so difficult to get gamma3 into simpler form, but it is the same as gamma2

simplify(gamma3 - gamma2)

ans = 

Get the numerical representation assuming a sampling period of T = 0.1

vpa(subs(phi,T,0.1),4)

ans = 

vpa(subs(gamma1,T,0.1),4)

ans = 

Show that phi and gamma can be computed numerically.

First approach

T = 0.1;
phi = expm(A*T)

phi = 2×2

0.9952 0.0950 -0.0950 0.9002

gamma1 = integral(@(tau) (expm(A*tau)*B),0,T,'ArrayValued',true)

gamma1 = 2×1

0.0048 0.0950

gamma2 = A\(phi - eye(2))*B

gamma2 = 2×1

0.0048 0.0950

temp = expm([A*T B*T;zeros(1,3)]);
phi3 = temp(1:2,1:2)

phi3 = 2×2

0.9952 0.0950 -0.0950 0.9002

gamma3 = temp(1:2,3)

gamma3 = 2×1

0.0048 0.0950

Answer 2

[Sam Chak]

@Rajesh Kanesan

Thanks for showing your calculation of the Gamma or Γ, I see now... Given the matrices and , you want to go from the continuous-time

to the discrete-time

where , to obtain

in terms of the sampling period T.

Since , and

if what I think about what you want is correct, then you may use this:

syms T
A = [0 1; -1 -1]
B = [0; 1]
Ad = expm(A*T)
Bd = A\(Ad - eye(size(A)))*B
Ad = simplify(Ad)
Bd = simplify(Bd)

You can check with @Paul if this approach can be admitted or not.

[Sam Chak]

@Rajesh Kanesan

If the sampling time T is given as shown in your assignment, you can verify the answer:

T = 0.1;
A = [0 1; -1 -1]
B = [0; 1]
Ad = expm(A*T)
Bd = A\(Ad - eye(size(A)))*B

Ad =

0.995167 0.095004

-0.095004 0.900163

Bd =

4.8334e-03

9.5004e-02

[Paul]

That's a cool equation for Bd. Have not seen that before.

[Rajesh Kanesan]

Thanks, @Sam Chak, I have submitted my assignment based on my manual calculation where I followed the steps that I was taught and tried to follow the same steps to answer another question but wished to validate my answers in MATLAB. My calculated phi value for the workout posted earlier is as shown below. step 2 was to find a Gamma matrix with the variable T and simulate the system with variable sampling intervals.

I replaced all the T with q as shown below. Could any one of you help, @Paul me to verify the answer with Matlab?

[Sam Chak]

Hi @Paul,

Thanks for showing the symbolic solution, which can be solved in a few ways. Actually at first I felt lost when doing this problem because my integration skills have rusted for a long time. Then I took it slowly.

[Paul]

At the risk of stating the obvious, need to enusre that A is invertible before using this equation in Matlab.

Also, something changed between 2022a and 2021b.

Running here on Answers with 2022a with A singular:

A = [1 1;0 0];
syms T
A\(expm(A*T-eye(2)))

Warning: Solution does not exist because the system is inconsistent.

ans = 

But when I run this same problem on my local installation of 2021b I get

>> A = [1 1; 0 0];

>> syms T

>> A\(expm(A*T) - eye(2))

Warning: Solution is not unique because the system is rank-deficient.

> In symengine

In sym/privBinaryOp (line 1136)

In \ (line 497)

ans =

[exp(T) - 1, exp(T) - 1]

[ 0, 0]

[Sam Chak]

Thanks @Paul for highlighting it. The formula only works if the matrix is non-singular.