matlab最小二乘辨识多项式系数.rar

clc clear close all % x—输入，y—输出，y_e—估计模型输出 % 估计模型为：y_e = a0 + a1*x + + 2*x^2 + a3*x^3 + a4*x^4 + a5*x^5 % 最小二乘估计目标函数： min 0.5 * Σ(y - y_e)^2 %% 读取Excel数据 data = xlsread('1000.xlsx','sheet1','A1:B5000'); N = 1000; x = data(:,1); % 输入量 y = data(:,2); % 输出量 n = 1000; %% 方法1：利用polyfit函数求解 % 构造计算5阶多项式的系数 p = polyfit(x,y,5); % 利用polyfit得到5阶多项式的系数 % 根据所求系数，计算多项式曲线的估计值 y1 = polyval(p,x); %% 方法2：构造偏导矩阵方程进行求解 % 对a0求偏导 H11 = N; H12 = sum(x); H13 = sum(x.^2); H14 = sum(x.^3); H15 = sum(x.^4); H16 = sum(x.^5); % 对a1求偏导 H21 = sum(x); H22 = sum(x.^2); H23 = sum(x.^3); H24 = sum(x.^4); H25 = sum(x.^5); H26 = sum(x.^6); % 对a2求偏导 H31 = sum(x.^2); H32 = sum(x.^3); H33 = sum(x.^4); H34 = sum(x.^5); H35 = sum(x.^6); H36 = sum(x.^7); % 对a3求偏导 H41 = sum(x.^3); H42 = sum(x.^4); H43 = sum(x.^5); H44 = sum(x.^6); H45 = sum(x.^7); H46 = sum(x.^8); % 对a4求偏导 H51 = sum(x.^4); H52 = sum(x.^5); H53 = sum(x.^6); H54 = sum(x.^7); H55 = sum(x.^8); H56 = sum(x.^9); % 对a5求偏导 H61 = sum(x.^5); H62 = sum(x.^6); H63 = sum(x.^7); H64 = sum(x.^8); H65 = sum(x.^9); H66 = sum(x.^10); % 构造H矩阵 H = [H11 H12 H13 H14 H15 H16; H21 H22 H23 H24 H25 H26; H31 H32 H33 H34 H35 H36; H41 H42 H43 H44 H45 H46; H51 H52 H53 H54 H55 H56; H61 H62 H63 H64 H65 H66;]; %6*6 % 构造Q矩阵 Q1 = sum(y); Q2 = sum(y .* x); Q3 = sum(y .* x.^2); Q4 = sum(y .* x.^3); Q5 = sum(y .* x.^4); Q6 = sum(y .* x.^5); Q = [Q1 Q2 Q3 Q4 Q5 Q6]'; % 计算A矩阵 A = H \ Q;% 相当于H逆*Q % 根据所求系数，计算多项式曲线的估计值 y2 =A(1)*y+ A(2)*x + A(3)*x.^2 + A(4)*x.^3 + A(5)*x.^4 + A(6)*x.^5; %% 计算样本方差 error1 = sum((y1 - y).^2)/(n-1); error2 = sum((y2 - y).^2)/(n-1); %% 在窗口显示 % 从低阶到高阶系数分别为： disp(strcat('polyfit 从高阶到低阶系数分别为:' , num2str(p'))); disp(strcat('最小二乘法 从低阶到高阶系数分别为:' , num2str(A'))); % 误差 disp(strcat('样本方差为:' , num2str(error1))); disp(strcat('样本方差为:' , num2str(error2))); %% 画图 figure(1) plot(y,'b'); hold on plot(y1,'r'); hold on plot(y2,'g'); legend('output','polyfit','最小二乘法','FontSize',10,'location','northeast');