Answer To: SYSC5807X_final_2023w is the assignment
Aditi answered on Apr 17 2023
Advance Mathematics
1.
% Define the problem parameters
c = [c1, c2, ..., cn]; % coefficient vector
A = eye(n); % identity matrix
b = [b1, b2, ..., bn]; % upper bound vector
lb = [a1, a2, ..., an]; % lower bound vector
% Solve the linear programming problem using linprog function
[x, fval] = linprog(c, A, b, [], [], lb);
% Display the optimal solution and objective function value
fprintf('Optimal solution: x* = [%s]\n', num2str(x));
fprintf('Objective function value: f(x*) = %f\n', fval);
2.
syms x1 x2 x3 lambda mu
% Define the objective function and constraints
f = [c1; c2; c3]' * [x1; x2; x3];
g = [a1 <= x1 <= b1, a2 <= x2 <= b2, a3 <= x3 <= b3, norm([x1;x2;x3]) <= 1];
% Write the Lagrangian function
L = f + lambda' * g;
% Find the gradient of the Lagrangian w.r.t. x and lambda
gradL_x = gradient(L, [x1, x2, x3]);
gradL_lambda = gradient(L, lambda);
% Write the KKT conditions
KKT_conditions = [gradL_x(1) == 0, gradL_x(2) == 0, gradL_x(3) == 0, gradL_lambda'*g
== 0, lambda >= 0, g >= 0, lambda'*g == 0];
Matlab code for (b) showing that if c - ATλ ≠ 0 where λ satisfies AATλ = Ac, then the
minimizer is given by x∗ = - (c - ATλ) / ||c - ATλ||:
% Check if c - AT*lambda != 0 where A*A' * lambda = A*c
if norm(c - A'*lambda) ~= 0
% Compute the minimizer
x_star = - (c - A'*lambda) / norm(c - A'*lambda);
else
disp('No minimizer exists since c - AT*lambda = 0')
end
3.
a.
Minimize ||Ax - b||∞:
% primal LP
cvx_begin
variable x(n)
variable t
minimize t
subject to
A*x - b <= t*ones(m,1)
A*x - b >= -t*ones(m,1)
cvx_end
b.
Minimize ||Ax - b||1 subject to ||x||∞ <= 1:
% primal LP
cvx_begin
variable x(n)
variable t
minimize ones(1,n)*abs(x-A\b)
subject to
max(abs(x)) <= t
cvx_end
c.
Minimize ||Ax - b||1 + ||x||∞:
% primal LP
cvx_begin
variable x(n)
variable t
minimize ones(1,n)*abs(x-A\b) + t
subject to
max(abs(x)) <= t
cvx_end
4.
a.
% Define problem data A = [1, 2; 3, 4]; % constraint matrix b = [5; 6]; % constraint
vector
% Define Lagrange dual function g = @(lambda) -1/4 * lambda' * A * A' * lambda +
lambda' * b;
% Test the dual function with some random values of lambda lambda = rand(2, 1); %
random vector of Lagrange multipliers g_lambda = g(lambda);
disp(['The value of the Lagrange dual function is ', num2str(g_lambda)])
b.
% Define problem data H = [2, -1; -1, 2]; % quadratic coefficient matrix p = [1; 2]; %
linear coefficient vector c = 3; % constant term P = [1, 0; 0, 1]; % quadratic constraint
matrix
% Define Lagrange dual function g = @(lambda) c - 1/4 * lambda' * (P * (H + lambda' *
P)^(-1) * P' + eye(2)) * lambda;
% Test the dual function with some random values of lambda lambda = rand(2, 1); %
random vector of Lagrange multipliers g_lambda = g(lambda);
disp(['The value of the Lagrange dual function is ', num2str(g_lambda)]);
5.
(a) Using forward Euler method with sampling rate 10 [Hz], we can discretize the
continuous-time state equation as:
x(k+1) = x(k) + dt * Ax(k) y(k) = Cx(k)
where dt = 1/10 is the time step, A = [1 1 -2; 2 0 0; 0 0 0] is the state matrix, and C = [1
0 0] is the output matrix.
To calculate the discretized system matrices Ad and Bd, we can use the following
formulas:
Ad = I + dt * A Bd = dt * B
where B = [0; p; 0] is the input matrix.
So, using the forward Euler method, we have:
Ad = eye(3) + dt * A Bd = dt * B
where B = [0; p; 0] is the input matrix.
To determine whether the discretized system is open-loop stable, we can check the
eigenvalues of Ad. If all the eigenvalues are inside the unit circle, then the system is
stable. We can use the MATLAB function eig(Ad) to find the eigenvalues.
(b) To design the controller for the discrete-time system, we first need to find the
discrete-time poles (eigenvalues) of the closed-loop system. We want the poles to be
at 2 and 1 +/- j1 to satisfy the design specifications. We can use the MATLAB function
c2d to discretize the controller design in the continuous-time domain to the discrete-
time domain.
Assuming a controller of the form u(k) = p * r(k) - K * x(k), we can use the MATLAB
function place to design the feedback gain K to place the desired poles. Then, we can
use the MATLAB function dcgain to find the feedforward gain p that ensures zero
steady-state error for a step input.
Here is the MATLAB code to implement the design and simulation of the discrete-time
system:
% Define problem data A = [1, 1, -2; 2, 0, 0; 0, 0, 0]; % state matrix B = [0; 1; 0]; %
input matrix C = [1, 0, 0]; % output matrix dt = 1/10; % sampling time
% Discretize the system using forward Euler method Ad = eye(3) + dt * A; Bd = dt * B;
% Check if the discretized system is open-loop stable eig_Ad = eig(Ad); if
all(abs(eig_Ad) < 1) disp('The discretized system is open-loop stable.'); else disp('The
discretized system is open-loop unstable.'); end
% Design the controller for the discrete-time system desired_poles = [2, 1+1j, 1-1j]; K
= place(Ad, Bd, desired_poles); p = -dcgain(Ad - BdK, BdC);
% Simulation of the closed-loop system with step inputs T = 0:dt:2; % simulation time
r = [5ones(1, length(T)/2), 2ones(1, length(T)/2+1)]; % step inputs x0 = [0; 0; 0]; %
initial state x = zeros(3, length(T)); % state trajectory y = zeros(1, length(T)); % output
trajectory u = zeros(1, length(T)); % input trajectory for k = 1:length(T) % Compute
control input u(k) = pr(k) - Kx(:,k);
% Simulate the system forward one step x(:,k+1) = Adx(:,k) + Bdu(k); y(k) = C*x(:,k);
end
% Plot the results subplot(2,1,1); plot(T, y); xlabel('Time (s)'); ylabel('Output y');
title('Step response of the closed-loop system');
subplot(2,1,2); plot(T, u); xlabel('Time (s)'); ylabel('Input u'); title('Control input of the
closed-loop system');
% Display the feedforward and feedback gains disp('The feedforward gain p is:');
disp(p); disp('The feedback gain K is:'); disp(K);
6.
(a) Force analysis:
The forces acting on the cart-pendulum system are the control force F, the
gravitational force mg, and the friction force bẋ. By applying Newton's second law and
torque balance to the cart and the pendulum, we get the following dynamic
equations:
Cart:
Mẍ = F - bẋ - mlsin(θ)θ̇²
Pendulum:
(I + ml²)θ̈ + bθ̇ + mgl sin(θ) = -m l cos(θ)ẍ
where M is the mass of the cart, m is the mass of the pendulum, b is the friction
coefficient, l is the length to the center of mass of the pendulum, I is the mass
moment of inertia of the pendulum, x is the position of the cart, θ is the angle of the
pendulum, F is the control force, g is the gravitational acceleration.
(b) State-space model:
Let x1 = x, x2 = ẋ, x3 = θ, and x4 = θ̇. Then we have the following state-space model:
ẋ = f(x, u) = Ax + Bu y = h(x) = Cx
where
A = [0 1 0 0; 0 -b/M -mg/M 0; 0 0 0 1; 0 -b/(Ml) -(M+m)g/(Ml) 0]
B = [0; 1/M; 0; 1/(M*l)]
C = [1 0 0 0; 0 0 1 0]
The matrices A and B are:
A = [0.0000 1.0000 0.0000 0.0000; 0.0000 -0.2000 -98.0000 0.0000; 0.0000 0.0000
0.0000 1.0000; 0.0000 -3.3333 -326.6667 0.0000]
B = [0.0000; 2.0000; 0.0000; 6.6667]
The system is controllable and observable since its controllability and observability
matrices have full rank:
rank([B AB A^2B A^3B]) = 4 rank([C; CA; CA^2; CA^3]) = 4
c.
Discretization and LQ-optimal control design:
The given continuous-time state-space model can be discretized using the exact
sampling method with sampling period T=0.1s. This can be done using the c2d
function in MATLAB. Then, a LQ-optimal controller can be designed using the dlqr
function.
Here is the complete MATLAB code for this part:
% Cart-pendulum system parameters
M = 0.5; m = 0.2; b = 0.1; l = 0.3; I = 0.006; g = 9.8;
% Continuous-time state-space model
Ac = [0 1 0 0; 0 -b/M g*m/M 0; 0 0 0 1; 0 b/(M*l) g*(M+m)/(M*l) 0];
Bc = [0; 1/M; 0; -1/(M*l)];
Cc = [1 0 0 0; 0 0 1 0];
Dc = [0; 0];
sysc = ss(Ac, Bc, Cc, Dc);
% Discrete-time state-space model
T = 0.1;
[Ad, Bd, Cd, Dd] = c2dm(Ac, Bc, Cc, Dc, T, 'zoh');
sysd = ss(Ad, Bd, Cd, Dd, T);
% LQ-optimal control design
Q = diag([10 1 100 1]);
R = 0.01;
[Klq, ~, ~] = dlqr(Ad, Bd, Q, R);
% Simulation time and initial conditions
t_sim = 5;
x0 = [1; 0; 0.3491; 0];
% State-feedback controller
K = Klq;
% Simulation with LQ-optimal controller
sim('cart_pendulum_simulink_LQ', t_sim);
% Plot results
figure;
subplot(2,1,1);
plot(tout, y(:,1), 'b-', 'LineWidth', 1.5);
hold on;...