%% 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   );
%       &uarr;  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



