Skip to content

Commit ea2494d

Browse files
MRGilakMRGilak
committed
added the DMC controller files
1 parent 02834ee commit ea2494d

26 files changed

Lines changed: 1093 additions & 6 deletions

‎DMC/compare_Ms.m‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
clear; close all; clc;
2+
3+
%% Set parameters
4+
params = get_base_params();
5+
M_values = [1, 6, 11];
6+
7+
%% Simulate
8+
for i = 1:length(M_values)
9+
params = get_base_params();
10+
params.M = M_values(i);
11+
params.R = 0.3107 * ones(params.M, 1);
12+
13+
% Simulate
14+
[results_lin{i}, sr_lin] = simulate_linear_system(params);
15+
[results_nonlin{i}, ~] = simulate_nonlinear_system(params);
16+
end
17+
18+
%% Plot results
19+
20+
plot_comparison_results(results_lin, results_nonlin, params, M_values, 'M', true)

‎DMC/compare_Ns.m‎

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,19 @@
1+
clear; close all; clc;
2+
3+
%% Set parameters
4+
params = get_base_params();
5+
N_values = [15, 55, 300];
6+
7+
%% Simulate
8+
for i = 1:length(N_values)
9+
params = get_base_params();
10+
params.N_model = N_values(i);
11+
12+
% Simulate
13+
[results_lin{i}, sr_lin] = simulate_linear_system(params);
14+
[results_nonlin{i}, ~] = simulate_nonlinear_system(params);
15+
end
16+
17+
%% Plot results
18+
19+
plot_comparison_results(results_lin, results_nonlin, params, N_values, 'N', true)

‎DMC/compare_Ps.m‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
clear; close all; clc;
2+
3+
%% Set parameters
4+
params = get_base_params();
5+
P_values = [11, 22, 44];
6+
7+
%% Simulate
8+
for i = 1:length(P_values)
9+
params = get_base_params();
10+
params.P = P_values(i);
11+
params.Q = ones(params.P, 1);
12+
13+
% Simulate
14+
[results_lin{i}, sr_lin] = simulate_linear_system(params);
15+
[results_nonlin{i}, ~] = simulate_nonlinear_system(params);
16+
end
17+
18+
%% Plot results
19+
20+
plot_comparison_results(results_lin, results_nonlin, params, P_values, 'P', true)

‎DMC/compare_Qs.m‎

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,24 @@
1+
clear; close all; clc;
2+
3+
%% Set parameters
4+
params = get_base_params();
5+
Q_values = [1, 10, 'decreasing'];
6+
7+
%% Simulate
8+
for i = 1:length(Q_values)
9+
params = get_base_params();
10+
if i < 3
11+
params.Q = Q_values(i) * ones(params.P, 1);
12+
end
13+
if i == 3
14+
params.Q = linspace(5, 0.01, params.P);
15+
end
16+
17+
% Simulate
18+
[results_lin{i}, sr_lin] = simulate_linear_system(params);
19+
[results_nonlin{i}, ~] = simulate_nonlinear_system(params);
20+
end
21+
22+
%% Plot results
23+
24+
plot_comparison_results_for_Q(results_lin, results_nonlin, params, Q_values, 'Q', true)

‎DMC/compare_Rs.m‎

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,19 @@
1+
clear; close all; clc;
2+
3+
%% Set parameters
4+
params = get_base_params();
5+
gamma_values = [0.01, 1, 100];
6+
7+
%% Simulate
8+
for i = 1:length(gamma_values)
9+
params = get_base_params();
10+
params.R = gamma_values(i) * 0.3107 * ones(params.M, 1);
11+
12+
% Simulate
13+
[results_lin{i}, sr_lin] = simulate_linear_system(params);
14+
[results_nonlin{i}, ~] = simulate_nonlinear_system(params);
15+
end
16+
17+
%% Plot results
18+
19+
plot_comparison_results(results_lin, results_nonlin, params, gamma_values, '\gamma', true)

‎DMC/compare_alphas.m‎

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,19 @@
1+
clear; close all; clc;
2+
3+
%% Set parameters
4+
params = get_base_params();
5+
alpha_values = [0.1, 0.5, 0.9];
6+
7+
%% Simulate
8+
for i = 1:length(alpha_values)
9+
params = get_base_params();
10+
params.alpha = alpha_values(i);
11+
12+
% Simulate
13+
[results_lin{i}, sr_lin] = simulate_linear_system(params);
14+
[results_nonlin{i}, ~] = simulate_nonlinear_system(params);
15+
end
16+
17+
%% Plot results
18+
19+
plot_comparison_results_for_alpha(results_lin, results_nonlin, params, alpha_values, true)

‎DMC/dmc_linear.m‎

Lines changed: 77 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,77 @@
1+
function Struct = dmc_linear(Struct)
2+
if ~isfield(Struct, 'upast') || isempty(Struct.upast)
3+
N = numel(Struct.sr);
4+
n = N - Struct.p;
5+
Struct.upast = zeros(n, 1);
6+
7+
% Ypast matrix
8+
step_response = Struct.sr(1:n);
9+
Struct.F = hankel(Struct.sr(2:Struct.p+1), Struct.sr(Struct.p+1:N)) - repmat(step_response(:)', Struct.p, 1); % Hankel
10+
11+
% Toeplitz Matrix
12+
Struct.G = toeplitz(Struct.sr(1:Struct.p), Struct.sr(1)*eye(1,Struct.m)); % Toeplitz
13+
14+
% DMC gain
15+
Q = diag(Struct.Q);
16+
R = diag(Struct.R);
17+
K = (Struct.G'*Q*Struct.G + R) \ (Struct.G'*Q);
18+
Struct.k = K(1, :);
19+
20+
% Control variables
21+
Struct.u = Struct.u0;
22+
Struct.u_prev = Struct.u0;
23+
Struct.d = 0;
24+
Struct.control = zeros(Struct.m,1);
25+
Struct.block_counter = 0;
26+
27+
if isfield(Struct, 'x0')
28+
Struct.y_prev = Struct.x0(1);
29+
Struct.y = Struct.x0(1);
30+
else
31+
Struct.y_prev = 0;
32+
Struct.y = 0;
33+
end
34+
end
35+
36+
% Disturbance
37+
if isfield(Struct, 'is_open_loop') && Struct.is_open_loop
38+
Struct.d = 0;
39+
else
40+
Struct.d = Struct.y - Struct.y_prev;
41+
end
42+
43+
D = ones(Struct.p,1) * Struct.d;
44+
45+
Ypast = Struct.F*Struct.upast + Struct.y;
46+
47+
% Clip the reference
48+
if Struct.is_programmed
49+
number_reference = numel(Struct.r);
50+
if number_reference >= Struct.p
51+
ref = Struct.r(1:Struct.p);
52+
else
53+
ref = [Struct.r(:); Struct.r(end) * ones(Struct.p - number_reference, 1)];
54+
end
55+
else
56+
ref = Struct.r(1) * ones(Struct.p, 1); % Unprogrammed: assume reference stays constant
57+
end
58+
59+
Struct.w = filter([0 (1-Struct.a)], [1 -Struct.a], ref, Struct.y); % Filtered Reference
60+
61+
Q = diag(Struct.Q);
62+
R = diag(Struct.R);
63+
K = (Struct.G'*Q*Struct.G + R) \ (Struct.G'*Q);
64+
Struct.Ebar = Struct.w - Ypast - D;
65+
Struct.control = K * Struct.Ebar;
66+
67+
% Apply control input
68+
Struct.du = Struct.control(1);
69+
Struct.upast = [Struct.du; Struct.upast(1:end-1)];
70+
Struct.u = Struct.u_prev + Struct.du;
71+
% Predict next output
72+
Struct.y_prev = Ypast(1) + Struct.G(1, 1) * Struct.du;
73+
74+
% Shift and update
75+
Struct.control = [Struct.control(2:end); 0];
76+
Struct.u_prev = Struct.u;
77+
end

‎DMC/dmc_nonlinear.m‎

Lines changed: 82 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,82 @@
1+
function Struct = dmc_nonlinear(Struct, x_now, ts)
2+
if ~isfield(Struct, 'upast') || isempty(Struct.upast)
3+
N = numel(Struct.sr);
4+
n = N - Struct.p;
5+
Struct.upast = zeros(n, 1);
6+
7+
% Ypast matrix
8+
step_response = Struct.sr(1:n);
9+
10+
% Toeplitz Matrix
11+
Struct.G = toeplitz(Struct.sr(1:Struct.p), Struct.sr(1)*eye(1,Struct.m)); % Toeplitz
12+
13+
% DMC gain
14+
Q = diag(Struct.Q);
15+
R = diag(Struct.R);
16+
K = (Struct.G'*Q*Struct.G + R) \ (Struct.G'*Q);
17+
Struct.k = K(1, :);
18+
19+
% Control variables
20+
Struct.u = -0.897;
21+
Struct.u_prev = -0.897;
22+
Struct.d = 0;
23+
Struct.control = zeros(Struct.m, 1);
24+
Struct.block_counter = 0;
25+
26+
if isfield(Struct, 'x0')
27+
Struct.y_prev = Struct.x0(1);
28+
Struct.y = Struct.x0(1);
29+
else
30+
Struct.y_prev = 0;
31+
Struct.y = 0;
32+
end
33+
end
34+
35+
% Disturbance
36+
if isfield(Struct, 'is_open_loop') && Struct.is_open_loop
37+
Struct.d = 0;
38+
else
39+
Struct.d = Struct.y - Struct.y_prev;
40+
end
41+
42+
D = ones(Struct.p,1) * Struct.d;
43+
44+
Ypast = zeros(Struct.p, 1);
45+
u_current = Struct.u_prev;
46+
x = x_now;
47+
for k = 1:Struct.p
48+
x = update_nonlinear_state(x, u_current, ts);
49+
Ypast(k) = x(1);
50+
end
51+
52+
if Struct.is_programmed
53+
number_reference = numel(Struct.r);
54+
if number_reference >= Struct.p
55+
ref = Struct.r(1:Struct.p);
56+
else
57+
ref = [Struct.r(:); Struct.r(end) * ones(Struct.p - number_reference, 1)];
58+
end
59+
else
60+
ref = Struct.r(1) * ones(Struct.p, 1); % Unprogrammed: assume reference stays constant
61+
end
62+
63+
Struct.w = filter([0 (1-Struct.a)], [1 -Struct.a], ref, Struct.y); % Filtered Reference
64+
65+
Q = diag(Struct.Q);
66+
R = diag(Struct.R);
67+
K = (Struct.G'*Q*Struct.G + R) \ (Struct.G'*Q);
68+
Struct.Ebar = Struct.w - Ypast - D;
69+
Struct.control = K * Struct.Ebar;
70+
71+
% Apply control input
72+
Struct.du = Struct.control(1);
73+
Struct.upast = [Struct.du; Struct.upast(1:end-1)];
74+
Struct.u = Struct.u_prev + Struct.du;
75+
76+
% Predict next output
77+
Struct.y_prev = Ypast(1) + Struct.G(1, 1) * Struct.du;
78+
79+
% Shift and update
80+
Struct.control = [Struct.control(2:end); 0];
81+
Struct.u_prev = Struct.u;
82+
end

‎DMC/general.m‎

Lines changed: 56 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,56 @@
1+
clear; close all; clc;
2+
3+
%% Define parameters
4+
params = struct();
5+
params.tf = 85;
6+
params.Ts = 0.2;
7+
params.N = params.tf / params.Ts;
8+
params.u0_nonlinear = -0.897;
9+
params.u0_linear = -1.7863;
10+
11+
bias = -1;
12+
span = 0.4;
13+
14+
%% Define input
15+
levels = bias * ones(1, 10) + span * [0, 0.6, 1.8, -0.3, -1.5, -0.6, 0, 2.4, -2.1, 0];
16+
spans = 1 * [5, 5, 7, 7, 7, 7, 6, 20, 12, 8];
17+
18+
span_steps = round(spans ./ params.Ts);
19+
ref_signal = [];
20+
for i = 1:length(levels)
21+
ref_signal = [ref_signal; repmat(levels(i), span_steps(i), 1)];
22+
end
23+
params.Ref = [ref_signal; repmat(levels(end), max(0, params.N - length(ref_signal)), 1)];
24+
25+
%% Define parameters
26+
params.P = 11;
27+
params.M = 8;
28+
params.alpha = 0.5;
29+
params.Q = 1 * ones(params.P, 1);
30+
params.R = 0.3107 * ones(params.M, 1);
31+
params.N_model = 55;
32+
params.disturbance_amp = 0.2;
33+
params.noise_power = -40;
34+
params.disturbance_start_time = 55;
35+
params.x0 = [-1; 0];
36+
params.is_programmed = true;
37+
params.is_open_loop = false;
38+
39+
%% Run simulations
40+
[results_lin, sr_lin] = simulate_linear_system(params);
41+
[results_nonlin, ~] = simulate_nonlinear_system(params);
42+
43+
Hankel = results_lin.F;
44+
Toeplitz = results_lin.G;
45+
gainMatrix = results_lin.gainMatrix;
46+
sr_all = results_lin.sr_all;
47+
48+
%% Plot results
49+
plot_simulation_results(results_lin, results_nonlin, params, true);
50+
% plot_step_response(params, sr_all);
51+
% plot_static_gain(-1);
52+
53+
54+
55+
56+

‎DMC/get_base_params.m‎

Lines changed: 38 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,38 @@
1+
function params = get_base_params()
2+
params = struct();
3+
params.tf = 85;
4+
params.Ts = 0.2;
5+
params.N = params.tf / params.Ts;
6+
params.u0_nonlinear = -0.897;
7+
params.u0_linear = -1.7863;
8+
9+
params.colors = {'b', 'r', 'g', 'm', 'c', 'k'};
10+
11+
bias = -1;
12+
span = 0.4;
13+
14+
% Define input
15+
levels = bias * ones(1, 10) + span * [0, 0.6, 1.8, -0.3, -1.5, -0.6, 0, 2.4, -2.1, 0];
16+
spans = 1 * [5, 5, 7, 7, 7, 7, 6, 20, 12, 8];
17+
18+
span_steps = round(spans ./ params.Ts);
19+
ref_signal = [];
20+
for i = 1:length(levels)
21+
ref_signal = [ref_signal; repmat(levels(i), span_steps(i), 1)];
22+
end
23+
params.Ref = [ref_signal; repmat(levels(end), max(0, params.N - length(ref_signal)), 1)];
24+
25+
% Define parameters
26+
params.P = 11;
27+
params.M = 8;
28+
params.alpha = 0.5;
29+
params.Q = 1 * ones(params.P, 1);
30+
params.R = 0.3107 * ones(params.M, 1);
31+
params.N_model = 55;
32+
params.disturbance_amp = 0.2;
33+
params.noise_power = -40;
34+
params.disturbance_start_time = 55;
35+
params.x0 = [-1; 0];
36+
params.is_programmed = true;
37+
params.is_open_loop = false;
38+
end

0 commit comments

Comments
 (0)