MFAPC代码（推导请看上几个博客）

链接上一个对该算法论文的复现，这里把相应的代码附上。

close all
clear all
N = 1000;
Ni = 4;
Ny = 8;
Nu = 4;
Np = 4;

a1  = 20;
b2 = 10;
b1 = 0.01; 
M = 1;

x_n = 2;
u_n = 2;
%x1 = zeros(x_n,N);
x = zeros(x_n,x_n,N);
u = zeros(u_n,N);

init_miu   = 0.2;
init_alpha = 0.2;
init_beta  = 0.2;
%sim('aim')

% r_mab = [dmab,ddmab];
% x_1_aim = dmab';
%x_aim = zeros(2,1000) + 1.5;
for k = 1:N
    if k<= 0.3*N
        x_aim(1,k) = 0.5;
        x_aim(2,k) = 1.5;
    else 
        x_aim(1,k) = 0.6;
        x_aim(2,k) = 1;
    end
    
end

T = zeros(x_n,u_n,N);
theta = zeros(x_n*Np,x_n,N);
phi= zeros(x_n*Np,u_n,N);

for i = 1:Ni
    %x(:,i) = [0.1;0.1;0.2;0.3;0.1;0.1];
    x(:,:,i) = [0.5,0;0.5,0];
    %x1(:,i) = zeros(x_n,1);
    u(:,i) = [0,0]';
    %T(:,:,i) = [0,0,0;0,0,0;0,0,0;0,0,0;0,0,0;0,0,0];
    T(:,:,i) = zeros(x_n,u_n)+[0.5,0.02;0.01,0.5];
end
phi_1 = [T(:,:,Ni)'];
for i = 1 : Np - 1
    phi_1 = [phi_1,T(:,:,Ni-i)'];
end
for i = 1:Ni
    phi(:,:,i) = phi_1';
end


for k = 1 : N
    a(k) = 1 + 0.1*sin(2*pi*k/1500);
    b(k) = 1 + 0.1*cos(2*pi*k/1500);
end
for k = Ni : N-Ny
    detla_y_k = x(:,1,k) - x(:,1,k-1);
    detla_u_k_1 = u(:,k-1) - u(:,k-2);
    T(:,:,k) = T(:,:,k-1) + (0.8*(detla_y_k-T(:,:,k-1)*detla_u_k_1)*detla_u_k_1')/(1 + sum(sum(detla_u_k_1.^2)) );
    theta(:,:,k) = theta(:,:,k-1) + phi(:,:,k-1)*((T(:,:,k)'-phi(:,:,k-1)'*theta(:,:,k-1)))/(0.9+sum(sum(phi(:,:,k-1).^2)));
    for m = 1 : x_n*Np
         for n = 1: x_n
              if abs(theta(m,n,k)) >= M
                    theta(m,n,k) = theta(m,n,1);
              end
         end
    end
    
    for j = 1 : Nu - 1
        phi_k = [T(:,:,k+j-1)'];
        for i = 1 : Np - 1
            phi_k = [phi_k,T(:,:,k+j-1-i)'];
        end
        phi(:,:,k+j-1) = phi_k';
    end
    
    for i = 1 : Nu - 1
         T(:,:,k+i) = theta(:,:,k)'*phi(:,:,k+i-1);
    end
    
    
%     for i = 1 : Nu 
%         for m = 1 : x_n
%             for n = 1 :u_n
%                 if m == n
%                     if abs(T(m,n,k+i-1)) <b2 || abs(T(m,n,k+i-1)) > a1*b2
%                         T(m,n,k+i-1) = T(m,n,1);
%                     end
%                 else
%                     if abs(T(m,n,k+i-1)) >b1 
%                         T(m,n,k+i-1) = T(m,n,1);
%                     end
%                 end
%             end
%         end
%         
%     end
    
    
    
    A_z = zeros(x_n,u_n);
    for i = 1:Ny
        A_1 = T(:,:,k);
        
        for j = 1:Nu-1
            if j+1 <= i
                A_1 = [A_1,T(:,:,k+j)];
            else
                A_1 = [A_1,A_z];
            end
        end
        if i >1 
            A = [A;A_1];
        else
            A = A_1;
        end
    end
    A_data(:,:,k) = A;
    Y_k = x(:,1,k);
    for i = 2:Ny
        Y_k = [Y_k;x(:,1,k)];
    end
    Y_aim_k_1 = x_aim(:,k+1);
    for i = 2:Ny
        Y_aim_k_1 = [Y_aim_k_1;x_aim(:,k+i)];
    end
    del_U = A' *(Y_aim_k_1-Y_k)/(4.5+sum(sum(A.^2))) ;
    del_u(:,k) = [del_U(1);del_U(2)];
    u(:,k) = u(:,k-1) + del_u(:,k);
    
    
    x(1,1,k+1) = x(1,1,k)^2/(1+x(1,1,k)^2) + 0.3*x(1,2,k);
    x(1,2,k+1) = x(1,1,k)^2/(1+x(1,2,k)^2+x(2,1,k)^2+x(2,2,k)^2)+a(k)*u(1,k);
    x(2,1,k+1) = x(2,1,k)^2/(1+x(2,1,k)^2) + 0.2*x(2,2,k);
    x(2,2,k+1) = x(2,1,k)^2/(1+x(1,1,k)^2+x(1,2,k)^2+x(2,2,k)^2)+b(k)*u(2,k);
  
end

for i = 1:N
    y1(i) = x(1,1,i);
    y2(i) = x(2,1,i);
    x3(i) = x(1,2,i);
    x4(i) = x(2,2,i);
end
y1_aim = x_aim(1,:);
y2_aim = x_aim(2,:);

figure
subplot(2,1,1)
plot(y1_aim)
hold on
plot(y1)
%ylim([0,2]);
legend('Target','Actual')
ylabel('Data')
rmse = mean((y1-y1_aim).^2);
title(sprintf('y1 RMSE=%.2f',rmse));
subplot(2,1,2)
error = y1-y1_aim;
scatter(1:1:N,error,'.')
xlabel('Time');ylabel('Error');
title('y1-误差图');

figure
subplot(2,1,1)
plot(y2_aim)
hold on
plot(y2)
%ylim([0,2.5]);
legend('Target','Actual')
ylabel('Data')
rmse = mean((y2-y2_aim).^2);
title(sprintf('y2 RMSE=%.2f',rmse));
subplot(2,1,2)
error = y2-y2_aim;
scatter(1:1:N,error,'.')
xlabel('Time');ylabel('Error');
title('y2-误差图');
% figure
% subplot(2,1,1)
% plot(u1)
% subplot(2,1,2)
% plot(u1)
% figure
% subplot(6,1,1)
% plot(x1)
% subplot(6,1,2)
% plot(u(1,:))
% subplot(6,1,3)
% plot(x2)
% subplot(6,1,4)
% plot(u(2,:))
% subplot(6,1,5)
% plot(x3)
% subplot(6,1,6)
% plot(x4)
% 
% 
% 
% figure
% subplot(2,1,1)
% plot(del_u(1,:))
% subplot(2,1,2)
% plot(del_u(2,:))

输出是

写在最后：

自从上一篇博客就有不少同学找我要代码，当时对MFAPC算法的研究仅仅是根据别人随便写的项目的本子做的，后来项目进度被人推翻重来MFAPC算法也被放弃，按理来说那时候就应该将代码分享出来的。但是那一年我的心情极其沉重阴暗，好在前段时间项目结题，最后也并未采用这一算法，因此今天将代码分享出来，希望能为每一个想毕业的牛马加一点夜草（马无夜草不肥）。