ECE 692: Discrete Time Modeling of Power Electronics
Buck-Boost Steady-State Simulation Example
The tabs below present 5 different approaches to solving the steady-state efficiency of the buck-boost converter modeled in Lectures 5-6. Each uses a slightly different method to solve the steady-state operation, and compute the averaged input and output power of the converter.
lsim Solution with Steady-State speedup
This code is the simplest approach to the problem, while still improving on the full lsim() approach that requires simulating over the initial transient to obtain the steady-state operating waveforms.
In this code, the 100 k$\Omega$ resistance $R_{out}$ is still present to overcome the singularity of ${\bf A}_1$. With every state matrix non-singular, the standard solution, ${\bf x}(t) = e^{{\bf A}t}{\bf x}_0 + {\bf A}^{-1}(e^{{\bf A}t}-{\bf I}){\bf B}{\bf u}$, still applies in every subinterval, so is used to solve for the steady-state initial condition ${\bf x}_0$. With this, lsim() only needs to cover one period to get the full steady-state waveform, rather than hundreds (or more) periods to wait for the initial transient to die out.
After using lsim() to get the full steady-state waveforms, trapz() is used to integrate the waveforms to find the averaged input and output power
tsteps = 1000; % Number of time-steps per period
%% Buck Boost Example
ron = 1e-3;
RL = 50e-3;
Co = .5e-6;
L = .3e-6;
fs = 1e6;
Ts = 1/fs;
D = 0.25;
Vg = 48;
Iload = 1;
Rout = 100e3; % Output resistor used to prevent A1 from being singular
K = [L, 0;
0, Co];
%% Define state space representations during each subinterval
% Subinterval I
A1p = [-(ron+RL), 0;
0, -1/Rout];
B1p = [1, 0;
0, -1];
A1 = K\A1p; % "\" computes the inverse of a matrix multiplied by another matrix/vector
B1 = K\B1p; % results in faster computation that inv(K)*A1p by not actually forming inv(K)
C1 = [ron, 0;
-ron, 1;
0, 1];
D1 = [0, 0;
1, 0;
0, 0];
% Subinterval II
A2p = [-(ron+RL), -1;
1, -1/Rout];
B2p = [0, 0;
0, -1];
A2 = K\A2p;
B2 = K\B2p;
C2 = [ron, 1;
-ron, 0;
0, 1];
D2 = [1, 0;
0, 0;
0, 0];
%% Steady-State Solution
u = [Vg; Iload];
ns = length(A1);
pi = length(u);
I = eye(ns);
eA1 = expm(A1*D*Ts); % state transition matrix during I
eA2 = expm(A2*(1-D)*Ts); % state transition matrix during II
X0 = (I-eA2*eA1)\( ...
eA2*(A1\(eA1-I))*B1*u + A2\(eA2-I)*B2*u );
% ↑ Note: parenthesis necessary to make sure the "\" doesnt take inverse of (eA2*A1)
XDTs = eA1*X0 + A1\(eA1-I)*B1*u; % Re-apply formula to find x(t=D*Ts)
XTs = eA2*XDTs + A2\(eA2-I)*B2*u; % Re-apply formula to find x(t=D*Ts)
%double-check to ensure steady-state was found.
assert(norm(XTs - X0) < 1e-9, 'Steady-State solution failed') %will throw an error if x0 =/= XTs
%% Reconstruct steady-state waveforms using lsim at solve X0
t = linspace(0, Ts, tsteps);
sys1 = ss(A1, B1, C1, D1); % create state space system during I
t1 = t(t<=D*Ts); % get timesteps during I
u1 = repmat(u, 1, length(t1)); % repeat u-vector (constant at each timestep)
[Y1,~,X1] = lsim(sys1,u1,t1,X0); % use lsim to get back time-domain waveform
sys2 = ss(A2, B2, C2, D2); % Repeat above for switching interval II
t2 = t(t>D*Ts);
u2 = repmat(u, 1, length(t2));
[Y2,~,X2] = lsim(sys2,u2,t2-t2(1),XDTs);
Y = [Y1; Y2]; % combine intervals I and II
X = [X1; X2];
%% Plot x(t) and y(t)
figure(1);
subplot(2,1,1)
plot(t*1e6, X, 'LineWidth', 3)
legend('{\it i_L} [A]', '{\it v_c} [V]');
ylabel('{\bf x}({\it{t}})')
subplot(2,1,2)
plot(t*1e6, Y, 'LineWidth', 3)
legend('{\it v_{ds1}} [V]', '{\it v_{ds2}} [V]', '{\it v_c} [V]');
xlabel('time [{\mu}sec]')
ylabel('{\bf y}({\it{t}})')
%% Calculate efficiency using lsim waveforms
ig = X(:,1).*(t<=D*Ts)';
vo = X(:,2);
Pin = trapz(t,Vg*ig)/Ts;
Pout = trapz(t,vo*Iload)/Ts;
eta = Pout/Pin
Time-Domain Solution to State Space
This code improves upon the lsim() example by using the direct solution to the state space equation within each interval to generat the waveforms for each of the states. As examined in HW3, this bypasses some of the limitations of time-stepping when generating waveforms, as the direct solution will remain correct at each time step, even if a very low sampling frequency is employed (i.e. a small value of tsteps).
However, with a small value of tsteps, the interpolation of the waveforms between the (correct-valued) sample points will be errant, so directly integrating the waveforms to calculation averaged powers will still require a large value of tsteps. This code additionally bypasses this issue by directly solving the integral of the state vector, rather than numerically integrating the waveform. As a result, even with a small sampling rate, the resultant calculation for $\eta$ will be correct
tsteps = 1000; % Number of time-steps per period
%% Buck Boost Example
ron = 1e-3;
RL = 50e-3;
Co = .5e-6;
L = .3e-6;
fs = 1e6;
Ts = 1/fs;
D = 0.25;
Vg = 48;
Iload = 1;
Rout = 100e3; % Output resistor used to prevent A1 from being singular
K = [L, 0;
0, Co];
%% Define state space representations during each subinterval
% Subinterval I
A1p = [-(ron+RL), 0;
0, -1/Rout];
B1p = [1, 0;
0, -1];
A1 = K\A1p; % "\" computes the inverse of a matrix multiplied by another matrix/vector
B1 = K\B1p; % results in faster computation that inv(K)*A1p by not actually forming inv(K)
C1 = [ron, 0;
-ron, 1;
0, 1];
D1 = [0, 0;
1, 0;
0, 0];
% Subinterval II
A2p = [-(ron+RL), -1;
1, -1/Rout];
B2p = [0, 0;
0, -1];
A2 = K\A2p;
B2 = K\B2p;
C2 = [ron, 1;
-ron, 0;
0, 1];
D2 = [1, 0;
0, 0;
0, 0];
%% Steady-State Solution
u = [Vg; Iload];
I = eye(ns);
eA1 = expm(A1*D*Ts); % state transition matrix during I
eA2 = expm(A2*(1-D)*Ts); % state transition matrix during II
X0 = (I-eA2*eA1)\( ...
eA2*(A1\(eA1-I))*B1*u + A2\(eA2-I)*B2*u );
% ? Note: parenthesis necessary to make sure the "\" doesnt take inverse of (eA2*A1)
XDTs = eA1*X0 + A1\(eA1-I)*B1*u; % Re-apply formula to find x(t=D*Ts)
XTs = eA2*XDTs + A2\(eA2-I)*B2*u; % Re-apply formula to find x(t=D*Ts)
%double-check to ensure steady-state was found.
assert(norm(XTs - X0) < 1e-9, 'Steady-State solution failed') %will throw an error if x0 =/= XTs
%% Reconstruct steady-state waveforms using direct solution
t = linspace(0, Ts, tsteps);
X = zeros(2,length(t)); %Pre-allocate
Y = zeros(3,length(t)); %Pre-allocate
for i = 1:1:length(t)
ti = t(i);
if(ti<= D*Ts)
X(:,i) = expm(A1*ti)*X0 + A1\(expm(A1*ti)-eye(2))*B1*u;
Y(:,i) = C1*X(:,i) + D1*u;
else
X(:,i) = expm(A2*(ti-D*Ts))*XDTs + A2\(expm(A2*(ti-D*Ts))-eye(2))*B2*u;
Y(:,i) = C2*X(:,i) + D2*u;
end
end
%% Plot x(t) and y(t)
figure(1);
subplot(2,1,1)
plot(t*1e6, X, 'LineWidth', 3)
legend('{\it i_L} [A]', '{\it v_c} [V]');
ylabel('{\bf x}({\it{t}})')
subplot(2,1,2)
plot(t*1e6, Y, 'LineWidth', 3)
legend('{\it v_{ds1}} [V]', '{\it v_{ds2}} [V]', '{\it v_c} [V]');
xlabel('time [{\mu}sec]')
ylabel('{\bf y}({\it{t}})')
%% Calculate efficiency using explicit integration of state space solution
intX1 = A1\((eA1-I)*X0 + A1\(eA1-I-A1*D*Ts)*B1*u); %Integral of states during t1
intX2 = A2\((eA2-I)*XDTs + A2\(eA2-I-A2*(1-D)*Ts)*B2*u); %Integral of states during t2
Vout = (1/Ts)*(intX1(2)+intX2(2)); %|Ts
Ig = (1/Ts)*intX1(1); %|Ts (just iL(t) in interval 1
Pout = Vout*Iload;
Pin = Vg*Ig;
eta = Pout/Pin
Augmented State Space (Single Term)
Though improved, the previous code still relied on an artificial $R_{out}$ resistance to enforce non-singular ${\bf A}_1$. In this code, and the following two, this resistor is removed and $|{\bf A}_1| = 0$.
The code first shows that the approach from the previous tabs fails in this case. All calculations still take place, but the first assertion shows that the resulting steady-state solution is not valid, because the calculated ${\bf x}_0$ is NaN.
Then this code uses the simplest implementation of the Augmented State Space approach. In this implementaiton, only the ${\bf A}_1^{-1}(e^{{\bf A}_1t_1}-{\bf I}){\bf B_1}{\bf u}$ term is calculated using the Augmented State Space, then the term is extracted from the resulting matrix and used in place of the incalculable term in the direct solution approach.
Because the direct integral of the state vector cannot be calculated when ${\bf A}_1$ is not invertible, this code goes back to the time-stepping approach, integrating the direct solution with trapz().
tsteps = 1000; % Number of time-steps per period
%% Buck Boost Example
ron = 1e-3;
RL = 50e-3;
Co = .5e-6;
L = .3e-6;
fs = 1e6;
Ts = 1/fs;
D = 0.25;
Vg = 48;
Iload = 1;
% Rout = 100e3; % No output resistor -> A1 is singular
K = [L, 0;
0, Co];
%% Define state space representations during each subinterval
% Subinterval I
A1p = [-(ron+RL), 0;
0, 0];
B1p = [1, 0;
0, -1];
A1 = K\A1p; % "\" computes the inverse of a matrix multiplied by another matrix/vector
B1 = K\B1p; % results in faster computation that inv(K)*A1p by not actually forming inv(K)
C1 = [ron, 0;
-ron, 1;
0, 1];
D1 = [0, 0;
1, 0;
0, 0];
% Subinterval II
A2p = [-(ron+RL), -1;
1, 0];
B2p = [0, 0;
0, -1];
A2 = K\A2p;
B2 = K\B2p;
C2 = [ron, 1;
-ron, 0;
0, 1];
D2 = [1, 0;
0, 0;
0, 0];
%% Steady-State Solution - Normal Approach
u = [Vg; Iload];
ns = length(A1);
pi = length(u);
I = eye(ns);
eA1 = expm(A1*D*Ts); % state transition matrix during I
eA2 = expm(A2*(1-D)*Ts); % state transition matrix during II
X0 = (I-eA2*eA1)\( ...
eA2*(A1\(eA1-I))*B1*u + A2\(eA2-I)*B2*u );
% ↑ Note: parenthesis necessary to make sure the "\" doesnt take inverse of (eA2*A1)
XDTs = eA1*X0 + A1\(eA1-I)*B1*u; % Re-apply formula to find x(t=D*Ts)
XTs = eA2*XDTs + A2\(eA2-I)*B2*u; % Re-apply formula to find x(t=D*Ts)
%double-check to ensure steady-state was found. Note: WILL fail in this example
try
assert(norm(XTs - X0) < 1e-9, 'Steady-State solution failed')
catch
display("Normal approach didn't work due to A1 being singular.");
display(X0);
end
%% Use Augmented State Space to calculate only the offending term
Aa1 = [A1, B1*u; zeros(1,ns+1)]; %Form Augmented matrix
eAa1 = expm(Aa1*D*Ts);
assert(norm(eAa1(1:2,1:2) - eA1) < 1e-9) %Double-check upper-left matrix is eA1
offendingTerm = eAa1(1:ns,ns+1); %Extract convolution integral term
% Replace term with singular A with Augmented State Space result.
X0 = (I-eA2*eA1)\( ...
eA2*offendingTerm + A2\(eA2-I)*B2*u );
XDTs = eA1*X0 + offendingTerm; % Re-apply formula to find x(t=D*Ts)
XTs = eA2*XDTs + A2\(eA2-I)*B2*u; % Re-apply formula to find x(t=D*Ts)
%Check solution (should pass here)
assert(norm(XTs - X0) < 1e-9, 'Augmented Steady-State solution failed');
%% Reconstruct steady-state waveforms using direct solution (With Aug SS)
t = linspace(0, Ts, tsteps);
X = zeros(2,length(t)); %Pre-allocate
Y = zeros(3,length(t)); %Pre-allocate
for i = 1:1:length(t)
ti = t(i);
if(ti<= D*Ts)
eAa1 = expm(Aa1*ti);
offendingTerm = eAa1(1:ns,ns+1); %Recalculate at the current time
X(:,i) = expm(A1*ti)*X0 + offendingTerm;
Y(:,i) = C1*X(:,i) + D1*u;
else
X(:,i) = expm(A2*(ti-D*Ts))*XDTs + A2\(expm(A2*(ti-D*Ts))-eye(2))*B2*u;
Y(:,i) = C2*X(:,i) + D2*u;
end
end
%% Plot x(t) and y(t)
figure(1);
subplot(2,1,1)
plot(t*1e6, X, 'LineWidth', 3)
legend('{\it i_L} [A]', '{\it v_c} [V]');
ylabel('{\bf x}({\it{t}})')
subplot(2,1,2)
plot(t*1e6, Y, 'LineWidth', 3)
legend('{\it v_{ds1}} [V]', '{\it v_{ds2}} [V]', '{\it v_c} [V]');
xlabel('time [{\mu}sec]')
ylabel('{\bf y}({\it{t}})')
%% Calculate efficiency using time-domain waveforms
ig = X(1,:).*(t<=D*Ts);
vo = X(2,:);
Pin = trapz(t,Vg*ig)/Ts;
Pout = trapz(t,vo*Iload)/Ts;
eta = Pout/Pin
Augmented State Space System Model #1
This code now uses the Augmented State Space method to alter the system being modeled. Instead of only calculating the single term which contains ${\bf A}_1^{-1}$, the system is modeled with a new homogeneous state-space description, $$\dot{\tilde{x}}(t) = \tilde{{\bf A}}_i{\tilde{\bf{x}}}(t)$$ where $$ \tilde{{\bf A}}_i = \left[ \begin{array}{c|c} {\bf A}_i & {\bf B}_i{\bf u}\\ \hline {\bf 0} & {\bf 0} \end{array} \right] $$ and $$ \tilde{{\bf x}}(t) = \left[ \begin{array}{c} {\bf x}(t) \\ \hline 1 \end{array} \right] $$ As discussed in Lecture 11, this description treats the inputs like an additional "state" which is constant within each subinterval. The result is a homogeneous system within each subinterval. The solution is therefore
$$\tilde{{\bf x}}(T_s) = (\prod_{i=k}^1e^{\tilde{{\bf A}}_it_i})\tilde{{\bf x}}_0$$Therefore, the steady-state initial condition $\tilde{{\bf x}}_0$ is in the nullspace of $\left( {\bf I}-\prod_{i=k}^1e^{\tilde{{\bf A}}_it_i}\right)$. The MATLAB command Z = null(M) returns a basis for the nullspace of matrix ${\bf M}$ (i.e. a vector ${\bf Z}$ such that ${\bf MZ}=0$). Of course, if ${\bf MZ}=0$, then ${\bf M}(k{\bf Z})=0$ for any scalar $k$, so we still need to find appropriate scaling for this basis such that the solution gives a valid $\tilde{{\bf x}}_0$. In this case, we know
$$ \tilde{{\bf x}}_0 = \left[ \begin{array}{c} {\bf x}_0 \\ \hline 1 \end{array} \right] $$
so we need only scale the nullspace basis such that the bottom-row element is unity.
Additionally, this code uses a further modification of the system to calculate the integrals of each state over each subinterval, which are then used to calculate $P_{out}$, $P_{in}$, and $\eta$, eliminating the need for time-stepping and numerical integration in the previous example. In addition to the $n_s$ original states and the one "state" capturing the effect of the inputs, $ns$ additional "states" are added with each additional state's derivative equal to the original states. The system is now $$\dot{\check{\bf x}}(t) = \check{{\bf A}}_i{\check{\bf{x}}}(t)$$ where the construction of each v-hat parameter is as follows $$ \left[ \begin{array}{c} \dot{\bf x}(t) \\ \hline 0 \\ \hline {\bf x}(t)\end{array} \right] = \left[ \begin{array}{c|c|c} {\bf A}_i & {\bf B}_i{\bf u} & {\bf 0}\\ \hline {\bf 0} & {\bf 0} & {\bf 0} \\ \hline {\bf I} & {\bf 0} & {\bf 0}\end{array} \right] \left[ \begin{array}{c} {\bf x}(t) \\ \hline 1 \\ \hline {\bf 0} \end{array} \right]$$ within the $i$th subinterval, with solution $${\check{\bf x}}(t) = \left[ \begin{array}{c} {\bf x}(t) \\ \hline 1 \\ \hline \int_0^t{\bf x}(\tau)d\tau \end{array} \right] = e^{\check{{\bf A}}_it_i}\check{{\bf x}}_0 $$ Therefore, by simply taking the matrix exponential across both subintervals, the final term in ${\check{\bf x}}(t)$ gives the integral of each of the states.
%% Buck Boost Example
ron = 1e-3;
RL = 50e-3;
Co = .5e-6;
L = .3e-6;
fs = 1e6;
Ts = 1/fs;
D = 0.25;
Vg = 48;
Iload = 1;
% Rout = 100e3; % No output resistor -> A1 is singular
K = [L, 0;
0, Co];
%% Define state space representations during each subinterval
% Subinterval I
A1p = [-(ron+RL), 0;
0, 0];
B1p = [1, 0;
0, -1];
A1 = K\A1p; % "\" computes the inverse of a matrix multiplied by another matrix/vector
B1 = K\B1p; % results in faster computation that inv(K)*A1p by not actually forming inv(K)
C1 = [ron, 0;
-ron, 1;
0, 1];
D1 = [0, 0;
1, 0;
0, 0];
% Subinterval II
A2p = [-(ron+RL), -1;
1, 0];
B2p = [0, 0;
0, -1];
A2 = K\A2p;
B2 = K\B2p;
C2 = [ron, 1;
-ron, 0;
0, 1];
D2 = [1, 0;
0, 0;
0, 0];
%% Steady-State Solution - Normal Approach
u = [Vg; Iload];
ns = length(A1);
pi = length(u);
I = eye(ns);
eA1 = expm(A1*D*Ts); % state transition matrix during I
eA2 = expm(A2*(1-D)*Ts); % state transition matrix during II
X0 = (I-eA2*eA1)\( ...
eA2*(A1\(eA1-I))*B1*u + A2\(eA2-I)*B2*u );
% ↑ Note: parenthesis necessary to make sure the "\" doesnt take inverse of (eA2*A1)
XDTs = eA1*X0 + A1\(eA1-I)*B1*u; % Re-apply formula to find x(t=D*Ts)
XTs = eA2*XDTs + A2\(eA2-I)*B2*u; % Re-apply formula to find x(t=D*Ts)
%double-check to ensure steady-state was found. Note: WILL fail in this example
try
assert(norm(XTs - X0) < 1e-9, 'Steady-State solution failed')
catch
fprintf(['Normal approach did not work due to A1 being singular; ', ...
'Attempt returned X0 = [%2.2f %2.2f]\n'], X0(1), X0(2));
end
%% Augmented State Space #1
Aa1 = [A1, B1*u; zeros(1,ns+1)]; %Form Augmented matrices
Aa2 = [A2, B2*u; zeros(1,ns+1)];
x0b = null(eye(ns+1) - ... %Find nullspace basis
expm(Aa2*(1-D)*Ts)*expm(Aa1*D*Ts));
X0 = x0b/x0b(end); %Scale basis so that final entry is "1"
XDTs = expm(Aa1*D*Ts)*X0;
XTs = expm(Aa2*(1-D)*Ts)*XDTs;
%Check solution (should pass here)
assert(norm(XTs - X0) < 1e-9, 'Augmented Steady-State solution failed');
%% Calculate efficiency using Augmented State Integration
Ai1 = [Aa1, zeros(ns+1,ns); I, zeros(ns,ns+1)]; %Form Augmented matrices
Ai2 = [Aa2, zeros(ns+1,ns); I, zeros(ns,ns+1)]; %including integration "dummy" states
X0i = [X0; zeros(ns,1)]; %append zeros for initial value of integral
%Compute average of states over full period
XTsi = expm(Ai2*(1-D)*Ts)*expm(Ai1*D*Ts)*X0i;
avgX = XTsi(end+1-ns:end)/Ts;
%Compute average of states over just the first subinterval
XDTsi = expm(Ai1*D*Ts)*X0i;
avgX_DTS = XDTsi(end+1-ns:end)/Ts;
Vout = avgX(2);
IL = avgX(1);
Ig = avgX_DTS(1);
Pout = Vout*Iload;
Pin = Ig*Vg;
eta = Pout/Pin
Augmented State Space System Model #2
This example closely follows the previous, but uses the Augmented system defined by $$ \tilde{{\bf A}}_i = \left[ \begin{array}{c|c} {\bf A}_i & {\bf B}_i\\ \hline {\bf 0} & {\bf 0} \end{array} \right] $$ and $$ \tilde{{\bf x}}(t) = \left[ \begin{array}{c} {\bf x}(t) \\ \hline {\bf u} \end{array} \right] $$ Because the system has $p_i=2$ inputs, the calculation of the nullspace to solve for $ \tilde{{\bf x}}_0$ returns $p_i=2$ bases. As before, we need to find an appropriate linear combination of these bases to yield a valid solution. With the two bases as $ {\bf Z}_1$ and $ {\bf Z}_2$, the solution we want to find is $$ k_1{\bf Z}_1 + k_2{\bf Z}_2 = \tilde{{\bf x}}_0 = \left[ \begin{array}{c} {\bf x}_0 \\ \hline {\bf u} \end{array} \right]$$ with $k_1$ and $k_2$ scalar. We are trying to find {\bf x}_0, but we know {\bf u}, so we can construct a linear system from only the bottom $p_i$ rows of this equation $$ k_1\left[ \begin{array}{c} {\bf Z}_{1,1} \\ \vdots \\ {\bf Z}_{1,n_s} \\ \hline {\bf Z}_{1,n_s+1} \\ \vdots \\ {\bf Z}_{1,n_s+p_i} \end{array} \right] + k_2\left[ \begin{array}{c} {\bf Z}_{2,1} \\ \vdots \\ {\bf Z}_{2,n_s} \\ \hline {\bf Z}_{2,n_s+1} \\ \vdots \\ {\bf Z}_{2,n_s+p_i} \end{array} \right] = \tilde{{\bf x}}_0 = \left[ \begin{array}{c} {\bf x}_0 \\ \hline {\bf u} \end{array} \right]$$ These bottom rows constitute the equation, for the $p_i=2$ case $$ \left[ \begin{array}{c c} {\bf Z}_{1,n_s+1} & {\bf Z}_{2,n_s+1} \\ {\bf Z}_{1,n_s+p_i} & {\bf Z}_{2,n_s+p_i}\end{array} \right] \left[ \begin{array}{c} k_1 \\ k_2\end{array} \right] = {\bf u} $$ Therefore, $$\tilde{{\bf x}}_0 = \left[ \begin{array}{c | c} {\bf Z}_{1} & {\bf Z}_{2}\end{array} \right] \left[ \begin{array}{c c} {\bf Z}_{1,n_s+1} & {\bf Z}_{2,n_s+1} \\ {\bf Z}_{1,n_s+p_i} & {\bf Z}_{2,n_s+p_i}\end{array} \right]^{-1}{\bf u}$$
The remainder of the code uses the same approach to Augmented state integration as the previous example, with differences only due to the size of $\tilde{{\bf A}}_i$ in this case.
%% Buck Boost Example
ron = 1e-3;
RL = 50e-3;
Co = .5e-6;
L = .3e-6;
fs = 1e6;
Ts = 1/fs;
D = 0.25;
Vg = 48;
Iload = 1;
% Rout = 100e3; % No output resistor -> A1 is singular
K = [L, 0;
0, Co];
%% Define state space representations during each subinterval
% Subinterval I
A1p = [-(ron+RL), 0;
0, 0];
B1p = [1, 0;
0, -1];
A1 = K\A1p; % "\" computes the inverse of a matrix multiplied by another matrix/vector
B1 = K\B1p; % results in faster computation that inv(K)*A1p by not actually forming inv(K)
C1 = [ron, 0;
-ron, 1;
0, 1];
D1 = [0, 0;
1, 0;
0, 0];
% Subinterval II
A2p = [-(ron+RL), -1;
1, 0];
B2p = [0, 0;
0, -1];
A2 = K\A2p;
B2 = K\B2p;
C2 = [ron, 1;
-ron, 0;
0, 1];
D2 = [1, 0;
0, 0;
0, 0];
%% Steady-State Solution - Normal Approach
u = [Vg; Iload];
ns = length(A1);
pi = length(u);
I = eye(ns);
eA1 = expm(A1*D*Ts); % state transition matrix during I
eA2 = expm(A2*(1-D)*Ts); % state transition matrix during II
X0 = (I-eA2*eA1)\( ...
eA2*(A1\(eA1-I))*B1*u + A2\(eA2-I)*B2*u );
% ↑ Note: parenthesis necessary to make sure the "\" doesnt take inverse of (eA2*A1)
XDTs = eA1*X0 + A1\(eA1-I)*B1*u; % Re-apply formula to find x(t=D*Ts)
XTs = eA2*XDTs + A2\(eA2-I)*B2*u; % Re-apply formula to find x(t=D*Ts)
%double-check to ensure steady-state was found. Note: WILL fail in this example
try
assert(norm(XTs - X0) < 1e-9, 'Steady-State solution failed')
catch
fprintf(['Normal approach did not work due to A1 being singular; ', ...
'Attempt returned X0 = [%2.2f %2.2f]\n'], X0(1), X0(2));
end
%% Augmented State Space #2
Aa1 = [A1, B1; zeros(pi,ns+pi)]; %Form Augmented matrices
Aa2 = [A2, B2; zeros(pi,ns+pi)];
x0b = null(eye(ns+pi) - ... %Find nullspace bases (2, because length(u) = 2)
expm(Aa2*(1-D)*Ts)*expm(Aa1*D*Ts));
sf = x0b(ns+1:end,:)\u; %Solve linear combination that results
%in u occupying final entries of
%steady-state solution
X0 = x0b*sf;
XDTs = expm(Aa1*D*Ts)*X0;
XTs = expm(Aa2*(1-D)*Ts)*XDTs;
%Check solution (should pass here)
assert(norm(X0(ns+1:end) - u) < 1e-9, 'Augmented Steady-State solution failed');
assert(norm(XTs - X0) < 1e-9, 'Augmented Steady-State solution failed');
%% Calculate efficiency using Augmented State Integration
Ai1 = [Aa1, zeros(ns+pi,ns); I, zeros(ns,ns+pi)]; %Form Augmented matrices
Ai2 = [Aa2, zeros(ns+pi,ns); I, zeros(ns,ns+pi)]; %including integration "dummy" states
X0i = [X0; zeros(2,1)]; %append zeros for initial value of integral
%Compute average of states over full period
XTsi = expm(Ai2*(1-D)*Ts)*expm(Ai1*D*Ts)*X0i;
avgX = XTsi(end+1-ns:end)/Ts;
%Compute average of states over just the first subinterval
XDTsi = expm(Ai1*D*Ts)*X0i;
avgX_DTS = XDTsi(end+1-ns:end)/Ts;
Vout = avgX(2);
IL = avgX(1);
Ig = avgX_DTS(1);
Pout = Vout*Iload;
Pin = Ig*Vg;
eta = Pout/Pin