% Linierisasi
    % Model Ruang Keadaan dari Kendalian:
        % x1dot = x2
        % x2dot = (F - (0.5*m*g*sin(2*x3)) + (m*l*(x4^2)*sin(x3)))/(M + m*(sin(x3))^2)
        % x3dot = x4
        % x4dot = ((g*sin(x3)) - (x2dot*cos(x3)))/l
    % Yang harus dilinierisasi adalah x2dot dan x4dot
    % Linierisasi dilakukan dengan meng-assumsi-kan:
        % sudut teta sangat kecil ---> x3 mendekati 0 (posisi pendulum
            % tegak)
        % kecepatan sudut teta kecil ---> x4 mendekati 0 (pendulum hampir
            % tidak bergoyang)
    % Konsekuensi: (1) sin(x3) ----> x3
    %              (2) cos(x3) ----> 1
    %              (3) kuadrat angka kecil ---> 0
    % Jadi: 
    % x2dot = (F - (m*g*x3))/M
    %       = - (m*g/M)*x3 +(1/M)*F
    % x4dot = ((g*x3) - x2dot)/l
    %       = ((g*x3) - (-m*g/M)*x3 +(1/M)*F)/l
    %       = ((M*g +m*g)/Ml)*x3 - (1/(M*l))*F
    %       = ((M +m)*g/(M*l))*x3 - (1/(M*l))*F
    % Model Ruang Keadaan hasil linierisasi:
    %  xdot = Ax + Bu, u adalah F
g = 10 ; % percepatan gravitasi 10 m/sec2
l = 0.6 ; % panjang batang 60 cm
M = 0.25 ; % Massa kereta 250 gram
m = 0.025 ; % massa pendulum 25 gram
A = [0         1          0          0;
     0         0       -(m*g/M)      0;
     0         0          0          1;
     0         0      (M+m)*g/(M*l)  0]
B = [0; 0; 0; -1/(M*l)]
pause
% Kestabilan
lambda = eig(A)
pause
% Keterkendalian dan Keteramatan
% Matrix Keterkendalian:
P = [B A*B A*A*B A*A*A*B]
p = rank(P)
pause
C = [1 0 0 0; % output posisi kereta
     0 1 0 0; % output kecepatan kereta
     0 0 1 0; % output posisi sudut
     0 0 0 1] % output kecepatan sudut kereta
pause
% Matrix Keteramatan:
Q = [C' A'*C' A'*A'*C' A'*A'*A'*C']
q = rank(Q)
pause
% Penempatan Nilai Eigen
% menggunakan Umpan-balik Peubah Keadaan (State Variable Feedback)
% Nilai eigen matrik [A-B*K] harus ditempatkan di sebelah kiri 
    % sumbu khayal agar sistem stabil, misalnya
    % lambda1 = -10
    % lambda2 = -5
    % lambda3 dan 4 = -0.5 +/- i*4
lambdaBaru = [-10;-5;-0.5+(i*4);-0.5-(i*4)]
pause
K = place(A,B,lambdaBaru)
pause
% Pembulatan nilai K -----> Kb
Kb = [120 45 -15 -2.5]
pause
% Hitung ulang nilai-eigen matrix [A - B*Kb]
eig(A-B*Kb)
% Selesai
% IMPLEMENTASI nilai Gain Matrix Kb ke Kendalian aslinya