% ===========================================================
%====================  28/12/2023 ===========================
% ===========================================================

clc
clear all
close all
clear

%=========  Component Values and Input Parameters ===========
% ===========================================================
VO_ = 30; % Output voltage
D_ = 0.6; % Duty cycle
VG_ = 12; % Input voltage
C_ = 100e-6; % Capacitor
rC_ = 50e-3; % Capacitor ESR
L_ = 120e-6; % Indcutor
rL_ = 10e-3; % Indcutor ESR
R_ = 50; % Output load



syms  vg rg d rL L rC C R vC iL rds rD vD D IL VC VG VD s

%%%%%%%%%%%%%%%%%%%%  STEP_1 %%%%%%%%%%%%%%%%%%%%
% Identify the amount of capacitor and the inductor on the DC-DC converter circuit. 
% Each of these energy storage elements will be the space vector of the matrix.

%%%%%%%%%%%%%%%%%%%%  STEP_2 %%%%%%%%%%%%%%%%%%%%
% Write down the equation for each state vector for each half cycle.

%%%%%%%%%%%%%%%%%%%%  STEP_3 %%%%%%%%%%%%%%%%%%%%
% Transform the state equations to a compact matrix form.
% ___________during DTs______________.
% 
% d[x]/dt = A1*x(t) + B1*u(t)
% y1(t)    = C1*x(t) + E1*u(t)
% 
% ___________during (1-D)Ts______________
% 
% d[x]/dt = A2*x(t) + B2*u(t)
% y2(t)    = C2*x(t) + E2*u(t)
% 
% Matrix x(t) consist of the inductor current (iL) and capacitor voltage (vC)
% Matrix u(t) consist of input/controlable parameters 
% e.g. input voltage (vg),input current(ig), 
% Matrix y(t) consist of output/ control target parameters
% e.g. output voltage (vo) and output current (io).

%_________Transforming the equations to matrices___________.

A1 = [-(rL)/L 0;
    0 -1/(C*(R+C))];

A2 = [-(rL+(R*rC)/(R+rC))/L -(R/(R+rC))/L;
    (R/(R+rC))/C -(1/(R+rC))/C];

B1 = [1/L;
    0];

B2 = [1/L;
    0];

C1 = [0 R/(R+rC)];

C2 = [R*rC/(R+rC) R/(R+rC)];

E1 =[0];

E2 = [0];

X = [IL;
    VC];

U = [VG];

%%%%%%%%%%%%%%%%%%%%  STEP_4 %%%%%%%%%%%%%%%%%%%%
% Calculating the steady-state operating point
% It is done by multiplying the first half-cycle matrices (A1, B1, C1, E1) with D
% Then multiplying the second half-cycle matrices (A2, B2, C2, E2) with (1-D).
%
% ____note____
% The MATLAB function "simplify" will do some algebric simplifications after finish calculating.

A = simplify(D*A1 + (1-D)*A2);
B = simplify(D*B1 + (1-D)*B2);
CC =simplify(D*C1 + (1-D)*C2); 
E = simplify(D*E1 + (1-D)*E2);

%%%%%%%%%%%%%%%%%%%%  STEP_5 %%%%%%%%%%%%%%%%%%%%
% During steady-state the inductor current and the capacitor voltage are
% constant (diL/dt = 0, dvC/dt =0)
% So we can write the matrices to something like below
%  ________________________________________                
%  _____          0 = AX + BU          ____
%  _____          Y = CX +EU           ____ 
%  ________________________________________
% 
% Moving the AX to the left side, -AX = BU
% Isolating X so we can get X = -inv(A)BU
% Finaly, we substitute X to Y eaquation
% Y = C(-inv(A)BU) + EU or we write as Y = (-Cinv(A)B +E)U 
% ___________________________________________________                
%  _____          X = -inv(A)BU                  ____
%  _____          Y = C(-inv(A)BU) +EU           ____ 
%  __________________________________________________
% 

DC_Solution_X = -inv(A)*B*U ;
DC_Solution_Y = -CC*inv(A)*B*U + E*U ; % we name the matrix as CC in order to distinguish matrix C and capacitor C

disp('Operating point of converter')
disp('-------------------------')
disp('IL(A)=')
disp(simplify(DC_Solution_X(1,1)))
disp('VC(V)=')
disp(simplify(DC_Solution_X(2,1)))
disp('VO(V)=')
disp(simplify(DC_Solution_Y))
disp('VO_simplify(V)=')
disp(simplify(subs((DC_Solution_Y),[rL rC],[0 0])))
disp('----------------------------')

%%%%%%%%%%%%%%%%%%%%  STEP_6 %%%%%%%%%%%%%%%%%%%%
% Calculating the averaged small-signal state-space model
% The small signal model can be obtained by adding small perturbations to the averaged state-space model
% Please refere to the post on the website, 
% to know how the small perturbation calculations are done for this converter
% 
% ________Calculating the additional matrices (F and G) due to small perturbation on the duty cycle variable
% 
% !!!!!!!!!___The matrix Xx consist of DC solution for IL and VC.___!!!!!!!!

Xx = DC_Solution_X ;

F = A1*Xx + B1*U - (A2*Xx + B2*U);
G = C1*Xx + E1*U - (C2*Xx + E2*U);

% Obtaining the control-to-output (V_od) transfer function by setting the input voltage u(s) to 0
% Obtaining the input-to-output (V_oi) transfer function by setting the duty-cyle d(s) to 0
%  ______________________________________________________               
%  _____          V_od(s) = C inv(sI-A)F + G         ____
%  _____          V_oi(s) = C inv(sI-A)B + E         ____ 
%  ______________________________________________________

I = eye(2); % identity matrix

Vo_d = simplify((CC*(inv(s*I-A))*F + G),'Steps',50) % The transfer function equation in "VG, R, D, L, C, rC, rL, s"

Vo_d_ = simplify(subs((CC*(inv(s*I-A))*F + G),[VG D R L C rL rC],[VG_ D_ R_ L_ C_ rL_ rC_]),'Steps',50) 

%%%%%%%%%%%%%%%%%%%%  STEP_7 %%%%%%%%%%%%%%%%%%%%
% Plotting the transfer function into a Bode Plot

s=tf('s');
opts = bodeoptions ('cstprefs');
opts.FreqUnits = 'Hz';
Vo_d__ =- ((1665418330830000000000*s)/334250917 - 166625808533983600000000000000/334585167917)/(500501001*s^2 + 225101996750*s + 6685018340000000) - 25015020000/334585167917;
bode(Vo_d__)
grid on
grid minor

%======================= END ====================================
%================================================================.