基于广义逆矩阵理论的分段线性拟合方法2

[x,fval,exitflag,output]=particleswarm(fun,nvars)
[x,fval,exitflag,output]=particleswarm(fun,nvars,lb,ub)
[x,fval,exitflag,output]=particleswarm(fun,nvars,lb,ub,options)

fun

nvars

lb

ub

option

x

fval

exitflag

global Data_x Data_y
[Data_s]= particleswarm(@pwlf_objfunc,nvars,ib,ub,options)

function RMSE=pwlf_objfunc(Data_s)
global Data_x Data_y
%这里输入定义RMSE的代码
end

function RMSE= pwlf_objfunc(Data_s,Data_x,Data_y)
%这里输入定义RMSE的代码
end

func_obj=@(x) pwlf_objfunc(x,Data_x,Data_y)
[Data_s]= particleswarm(@func_obj,nvars,ib,ub,options)

clear,clc
%**************************************************************************
%程序用于分段线性拟合 
%输入参数为待拟合数据以及分段断点数量
%注意 待拟合数据点必须以列向量的形式输入
%**************************************************************************

%***************************待拟合数据及参数输入***************************
Data_x=linspace(0,5*pi,101)';                  %待拟合数据点x
Data_y=sin(Data_x)+0.2*sin(5*pi*Data_x);       %待拟合数据点y     
Split_pt_num=5;                                %分段断点数量
%**************************************************************************

%**************************分段线性拟合计算********************************
%创建目标函数匿名函数句柄用于优化计算
plot_status="false";
func_obj=@(x) pwlf_objfunc(x,Data_x,Data_y,plot_status);

%基于粒子群算法优化断点布置位置
lb=zeros(Split_pt_num,1)+min(Data_x);          %变量的下限
ub=zeros(Split_pt_num,1)+max(Data_x);          %变量的上限
options=optimoptions("particleswarm","PlotFcn","pswplotbestf");
[Data_s,fval,exitflag,output]=...
    particleswarm(func_obj,Split_pt_num,lb,ub,options);
%**************************************************************************

%************************绘制分段线性拟合结果******************************
%创建目标函数匿名函数句柄用于输出分段线性拟合结果
plot_status="true";
[RMSE]=pwlf_objfunc(Data_s,Data_x,Data_y,plot_status);
%**************************************************************************

function [RMSE]=pwlf_objfunc(Data_s,Data_x,Data_y,plot_status)
%**************************************************************************
%Data_s为拟合区段断点 以行向量形式指定
%Data_x和Data_y为待拟合数据点 以列向量形式指定
%**************************************************************************
%***************************待拟合数据及参数输入***************************
Data_num=size(Data_x,1);                        %待拟合数据点总数
%剔除重复断点并按升序排列
Data_s=unique(Data_s)';
%剔除位于拟合点边界的断点(只有可能为断点的最小值和最大值)
if Data_s(1)==Data_x(1)
    if Data_s(end)==Data_x(end)
        Interval_num=size(Data_s,1)-1;
        Data_s=Data_s(2:end-1);
    else
        Interval_num=size(Data_s,1);
        Data_s=Data_s(2:end);
    end
else
    if Data_s(end)==Data_x(end)
        Interval_num=size(Data_s,1);
        Data_s=Data_s(1:end-1);
    else
        Interval_num=size(Data_s,1)+1;
    end
end
%判断剔除操作完成后是否出现无断点的情况
if Interval_num==1
    RMSE=Inf;         %若无断点 假设目标函数取值为无穷大
    return;
else
    Data_s=[Data_s;max(Data_x)];
end
%**************************************************************************

%********************确定待拟合数据点所在区段的索引值**********************
Interval_index=zeros(Data_num,1);
Start_index=1;
for i=1:Data_num
    for j=Start_index:Interval_num
        %区段索引为1的特殊情况
        if j==1
            if Data_x(i)<=Data_s(j)
                Interval_index(i)=j;
                break;
            end
            continue;
        end
        %区段索引不为1的一般情况
        if (Data_x(i)>Data_s(j-1)&&Data_x(i)<=Data_s(j))
            Interval_index(i)=j;
            Start_index=j;
            break;
        end
    end
end
%**************************************************************************

%***************************确定等价关系矩阵A******************************
%初始化矩阵A
Matrix_A=zeros(Data_num,2*Interval_num);
%计算矩阵A的分量
for i=1:Data_num                   %行循环
    Matrix_A(i,2*(Interval_index(i)-1)+1)=Data_x(i);
    Matrix_A(i,2*(Interval_index(i)-1)+2)=1;
end
%**************************************************************************

%***************************生成连续性条件矩阵*****************************
%初始化存放不同区段的连续性条件矩阵的元胞数组
Matrix_C=cell(Interval_num-1,1);
%循环创建不同区段的连续性条件矩阵
for j=1:Interval_num-1
    Matrix_C{j,1}=[-1/Data_s(j) 1 1/Data_s(j)];
end
%**************************************************************************

%**************************组装转换矩阵分矩阵******************************
%初始化存放不同区段的转换矩阵分矩阵的元胞数组
Matrix_T=cell(Interval_num-1,1);
%考虑i=1的特殊情况[T1]
Matrix_T{1}=zeros(2*Interval_num,2*Interval_num-1);
Matrix_T{1}(1,1:3)=Matrix_C{1};               %填充连续性条件矩阵[C]
Matrix_T{1}(2:end,:)=eye(2*Interval_num-1);   %填充单位矩阵
%考虑i>1的一般情况[Ti]
for j=2:Interval_num-1
    Matrix_T{j}=zeros(2*Interval_num-j+1,2*Interval_num-j);
    Matrix_T{j}(1:j-1,1:j-1)=eye(j-1);          %填充上侧单位矩阵
    Matrix_T{j}(j,j:j+2)=Matrix_C{j};           %填充连续性条件矩阵
    Matrix_T{j}(j+1:end,j:end)=eye(2*Interval_num-2*j+1);   %填充下侧单位矩阵
end
%**************************************************************************

%****************************计算分段拟合直线系数**************************
%计算转换矩阵[Pai_T]
Matrix_Pai_T=Matrix_T{1};
for j=2:Interval_num-1
    Matrix_Pai_T=Matrix_Pai_T*Matrix_T{j};     %通过连乘分矩阵获得转换矩阵
end
%计算独立变量系数计算矩阵[E]
Matrix_E=Matrix_A*Matrix_Pai_T;
%计算独立变量系数计算矩阵[E]的广义逆矩阵
Matrix_E_inv=pinv(Matrix_E);
%计算分段样条曲线的系数矩阵
Matrix_Coeff=Matrix_Pai_T*Matrix_E_inv*Data_y;
%**************************************************************************

%************************计算分段线性拟合的均方根误差**********************
Fitting_y=zeros(Data_num,1);
for i=1:Data_num
    Fitting_y(i)=Matrix_Coeff(2*(Interval_index(i)-1)+1)*Data_x(i)+...
        Matrix_Coeff(2*(Interval_index(i)-1)+2);
end
RMSE=sqrt((sum((Data_y-Fitting_y).^2))/Data_num);
%**************************************************************************

%*************************判断是否需要绘制函数图像*************************
if plot_status=="true"
    %初始化用于绘制分段拟合直线的元胞数组
    Fitting_x=cell(Interval_num,1);
    Fitting_y=cell(Interval_num,1);
    for i=1:Interval_num
        if i==1
            Fitting_x{i}=linspace(min(Data_x),Data_s(i),51)';
        else
            Fitting_x{i}=linspace(Data_s(i-1),Data_s(i),51)';
        end
    end
    %计算分段样条拟合直线
    for i=1:Interval_num
        %基于分段样条拟合直线的系数进行计算
        Fitting_y{i}=Matrix_Coeff(2*(i-1)+1)*Fitting_x{i}+...
            Matrix_Coeff(2*(i-1)+2);
    end
    %绘制分段样条直线
    figure;
    scatter(Data_x,Data_y,70,"Filled","Black");hold on;
    for i=1:Interval_num
        plot(Fitting_x{i},Fitting_y{i},"LineWidth",2,"Color","red");
    end
    xlabel("x");ylabel("y");title("分段线性拟合结果");hold off;grid on;
    %输出分段样条直线系数
    fprintf("分段样条直线为:y=ax+b\n");
    for i=1:Interval_num
        fprintf("第%d条拟合直线系数为:",i);
        fprintf("a=%10.5f,",Matrix_Coeff(2*(i-1)+1));
        fprintf("b=%10.5f,",Matrix_Coeff(2*(i-1)+2));
        if i==1
            fprintf("取值范围为:[%10.5f,%10.5f]\n",min(Data_x),Data_s(i));
        else
            fprintf("取值范围为:[%10.5f,%10.5f]\n",Data_s(i-1),Data_s(i));
        end
    end
else
    return;
end
%**************************************************************************
end