least_square_paramete_solution.m

控制行业中重要的least square parameter solution

% Tasks:
% 1. Data pre-treatment using both methods discussed in the classroom;
% 2. Least squares parameter solution;
% 3. Plot both impulse response and step response;
% 4. Generate the FIR model's output with screw speed as the input and
%    compare this predicted output with the actual measured process output;
% 5. Compare the estimation results for the two data pre-treatment methods;
% 6. Validate your programs using the testing process;
% 7. Discuss your observations on the estimation results by varying the
%    number of terms (or model order N);


% load assignment2 data;
 load assignment2;

% Task 6 (different data can be loaded by adding and deleting 
%           '%' in front of the sentences)
% load noise free test data;
% load test_data_no_noise;u=u';y=y'; % refter to Figure 10a,10b
% load test data with noise;
% load test_data_noise;u=u';y=y'; % refter to Figure 11a,11b

% size of parameter matrix n;
clc;
t_sample=120;
t_inter=1;

% Task 7; (different t_sample can be loaded by adding and deleting 
%          '%' in front of the sentences))
% t_sample=50; % refter to Figure 12
% t_sample=150; % refter to Figure 13
n=t_sample/t_inter;

% set input signal matrix a which will be 2500*120 
a=u; % u is the input signal array, from u(1) to u(2500);
i=1;
j=1;
c=size(u); % c=(2500,1);
for i=1:n-1
    a=[a [zeros(j,1);a(1:(c(1)-i),1)]];
    j=j+1;
end

% Task 1:
% 1.1 Removing steady state values, method 1: minus the steady state value;
y_ss=mean(y); % to get the average value of y
u_ss=mean(u); % to get the average value of u
y_k=y-y_ss; % remove the steady state value
u_k=u-u_ss; % remove the steady state value
figure;
subplot(211);
plot(u_k); % refter to Figure 1
grid;
xlabel('k');
ylabel('u_k=u-u_ss');
subplot(212);
plot(y_k); % refter to Figure 1
grid;
xlabel('k');
ylabel('y_k=y-y_ss');

% 1.2 Removing steady state values, method 2: minus the steady state value;
c=size(u); % c=(2500,1);
Yop_k=y; % define Yop(k);
Yop_k1=[mean(y);y(1:(c(1)-1))]; % define Yop(k-1);
Uop_k=u; % define Uop(k);
Uop_k1=[mean(u);u(1:(c(1)-1))]; % define Uop(k-1);
Y_k=Yop_k-Yop_k1; % diff Y(k)=Yop(k)-Yop(k-1);
U_k=Uop_k-Uop_k1; % diff U(k)=Uop(k)-Uop(k-1);
figure;
subplot(211); 
plot(U_k); % refter to Figure 2
grid;
xlabel('k');
ylabel('difference U(k)=Uop(k)-Uop(k-1)');
subplot(212);
plot(Y_k); % refter to Figure 2
grid;
xlabel('k');
ylabel('difference Y(k)=Yop(k)-Yop(k-1)');

% Task 2:
% 2.1: obtain the best parameter vector b1 (120*1) using method 1;
% using method 1, set input signal matrix a1 which will be 2500*120;
a1=u_k; % a1 is the matrix of u from which the steady state value has been removed;
i=1;
j=1;
c=size(u_k); % c=(2500,1);
for i=1:n-1
    a1=[a1 [zeros(j,1);a1(1:(c(1)-i),1)]];
    j=j+1;
end
% to obtain the best parameter vector b1 (120*1) using method 1;
b1=inv(a1'*a1)*a1'*y_k; 


% 2.2: obtain the best parameter vector b2 (120*1) using method 2;
a2=U_k; % a2 is the data matrix of the input difference U_k
i=1;
j=1;
c=size(U_k); % c=(2500,1);
for i=1:n-1
    a2=[a2 [zeros(j,1);a2(1:(c(1)-i),1)]];
    j=j+1;
end
b2=inv(a2'*a2)*a2'*Y_k; % to obtain the best parameter vector b2 (120*1) using method 2

figure; % to plot b1 and b2
subplot(211);
plot(b1); % refter to Figure 3
grid;
xlabel('k');
ylabel('b1') 
title('the best parameter vector b1 (method 1)');
subplot(212);
plot(b2); % refter to Figure 3
grid;
xlabel('k');
ylabel('b2');
title('the best parameter vector b2 (method 2)')

% Task 3:
% 3.1: plot step response for both method 1 and method 2
% step response (method 1)
c=size(u);
d=ones(c(1),1); % define a signal step signal;
i=1;
j=1;
for i=1:n-1
    d=[d [zeros(j,1);d(1:(c(1)-i),1)]];
    j=j+1;
end
y_step1=d*b1; % y_step1 is the step response
figure;
subplot(211);
plot(y_step1(1:125)); % refter to Figure 4
grid;
xlabel('k');
ylabel('y step1')
title('step response with b1 using method 1');

% step response (method 2)
y_step2=d*b2; % y_step2 is the step response using method 2
subplot(212);
plot(y_step2(1:125)); % refter to Figure 4
grid;
xlabel('k');
ylabel('y step2')
title('step response with b2 using method 2');


% 3.2: plot impulse response for both method 1 and method 2
% impulse response (method 1)
e=[1;zeros(c(1)-1,1)]; % define a signal step signal;
i=1;
j=1;
c=size(u);
for i=1:n-1
    e=[e [zeros(j,1);e(1:(c(1)-i),1)]];
    j=j+1;
end
y_impulse1=e*b1; % y_impulse 1 is the impulse response
figure;
subplot(211);
plot(y_impulse1(1:125)); % refter to Figure 5
grid;
xlabel('k');
ylabel('y impulse1');
title('impulse response with parameter vector b1 using method 1');
%b1-y_impulse1(1:120); %b1 should be equal to y_impulse1(1:125)

% impulse response (method2)
y_impulse2=e*b2; % y_step1 is the impulse response
subplot(212);
plot(y_impulse2(1:125)); % refter to Figure 5
grid;
xlabel('k');
ylabel('y impulse2')
title('impulse response with parameter vector b2 using method 2');


% Task 4
% the comparison between the predicted output y_pre with the actual measured process output y
% 4.1: method 1
y_k_pre1=a1*b1;% to get the prediction of y_k using b1;
y_pre1=y_k_pre1+y_ss; % the predicted output y_pre1;
diff_y1=y-y_pre1; % the difference between the prediction and actual measurement;
figure;
subplot(211);
plot(y);
grid;
xlabel('k');
ylabel('y');
title('the actual measured process output');
subplot(212);
plot(y_pre1); % refter to Figure 6
grid;
xlabel('k');
ylabel('y pre1');
title('the predicted output using method 1');


% 4.2: method 2
Y_k_pre2=a2*b2;% to get the prediction of y_k using b2; 
y_pre2=zeros(c(1),1); 
y_pre2(1)=Yop_k(1); % to derive predicted output y_pre2, we assume the first number of y_pre2 y(1);
for m=2:c(1)
    y_pre2(m)=Y_k_pre2(m)+y_pre2(m-1); % derive the predicted output y_pre2
end
%y_pre2=Y_k_pre2+y_ss; % the predicted output y_pre2; do i plus y_ss or anything else?
diff_y2=y-y_pre2; % the difference between the prediction and actual measurement;
%plot(diff_y2(1:2500));
figure;
subplot(211);
plot(y);
grid;
xlabel('k');
ylabel('y');
title('the actual measured process output');
subplot(212);
plot(y_pre2); % refter to Figure 7
grid;
xlabel('k');
ylabel('y pre2');
title('the predicted output using method 2');

% Task 5:
% from above, we get the difference between prediction and actual measurement 
% using both method 1 and method 2. Now we compare them.
figure;
subplot(211);
plot(y_pre1); % refter to Figure 8
grid;
xlabel('k');
ylabel('y pre1');
title('the predicted output using method 1');
subplot(212);
plot(y_pre2); % refter to Figure 8
grid;
xlabel('k');
ylabel('y pre2');
title('the predicted output using method 2');

figure;
plot(diff_y1-diff_y2); % we can see from this figure that prediction is closer to actual 
                       % measurement in method 1 than in method 2;
                       % refter to Figure 9
grid;
xlabel('k');
ylabel('diff y1-diff y2');
title('comparison of the difference to actual measurement using both methods');