首页 数学基础

  区间逐维分析方法借鉴了代理模型的思想，不同之处在于其利用参考点 (Reference Point) 处输入参数向量每个维度上的超平面截取结构真实响应所形成的曲面，并基于勒让德正交多项式构建所截得曲线的多项式近似模型，进一步可以计算在对应维度上的最大值点与最小值点。

  区间逐维分析方法主要使用了高斯采样和勒让德正交多项式拟合两个过程，大致流程图可参考下图。

# 一、输入逐维高斯采样

  区间逐维分析方法的第一步就是根据输入的维度依次进行高斯采样。在对第ii 维的参数进行高斯采样时，需要固定其他维度的参数为中心值，之后对第ii 维的参数进行高斯采样。其中，高斯采样的区间为 [-1, 1]。高斯采样的 matlab 代码如下（输入为 s 时，会返回 s+1 个高斯采样点）：

function y= Gauss(n)
    % gauss-legendre 求积法点集的求解程序 论文里原代码
    if n==0
        y=0;
    else
        beta = 0.5./sqrt(1-(2*(1:n)).^(-2));% B = sqrt (X) 返回数组 X 的每个元素的平方根
        T = diag(beta,1)+diag(beta,-1);% 取 beta 中的元素为上方第 1 条对角线和下方第一条对角线
        [~,D]=eig(T);% eig 为求特征值和特征向量
        x=diag(D);% 取特征值
        [x,~]=sort(x);
        y=x;
    end
end

  完成高斯采样后，需要将得到的各个标准采样点等比例放缩到响应参数的区间，得到各个真实采样点。之后对每一个真实采样点进行映射运算得到每一个真实采样点对应的目标函数值。对应的 matlab 代码如下

%% 计算样本点对应的函数值
for i=1:m  % 对每个优化参数进行一次区间逐维分析
    pC=pC_All(i);  % 第 i 个参数的中值
    pR=pR_All(i);  % 第 i 个参数的半径
    x_gauss=Gauss(s);  % 式 3.82 获取高斯积分点 有 s+1 个
    p_gauss=pC+pR.*x_gauss;  % 式 3.83 + 式 3.84 获取第 i 个参数的真实采样点
    
    % 式 3.85 + 式 3.86
    S_p_gauss=repmat(pC_All, s+1, 1);
    for j=1:s+1
        S_p_gauss(j,i)=p_gauss(j);
    end
    
    % 式 3.87 + 式 3.88
    for j=1:s+1
        u_all=MyFun(S_p_gauss(j,:));% 此时，解中包含其他维度的值
        cd(j+(s+1)*(i-1),:)=u_all;% 保存所有真实计算的结果 列是响应向量 行是
    end
end

# 二、勒让德多项式拟合

  针对第一章中采样得到的标准采样点，对每一个目标函数维度和输入参数维度，依次使用勒让德多项式进行拟合。勒让德多项式拟合的 matlab 代码如下：

u=cd((s+1)*(i-1)+1:(s+1)*i,q);% 对应参数和对应目标函数的函数关系
% 求解 n 阶 Legendre 拟合多项式
syms x
for k=0:n
    % 求 A_Tq_k_1
    A_Tq_k_1=0;
    for j=1:s+1
        if k==0
            L_k_x=1;
        else
            L_k(x)=Legendre_k(k,x);
            L_k_x=L_k(x_gauss(j));
        end
        L_s_diff(x)=diff(Legendre_k(s+1,x),1);
        A_Tq_k_1=A_Tq_k_1+((2*k+1)*u(j)*L_k_x)/((1-x_gauss(j)^2)*vpa((L_s_diff(x_gauss(j)))^2));% Lk (x) 的系数 公式 3.90
    end
    % 利用 k 和 A_Tq_k_1 求解 T (x)（每一次都累加，共累加 n 次，n 阶平方逼近）
    if k==0
        T(x)=A_Tq_k_1*Legendre_k(k,x);
    else
        T(x)=T(x)+A_Tq_k_1*Legendre_k(k,x);
    end
end

  其中 Legendre_k 表示勒让德多项式函数， k 是拟合阶数， x 代表符号变量（可以替代为 ~ ）。 Legendre_k 函数可见如下代码：

function L_k=Legendre_k(k,x)
    syms x
    L_k=1./double((2.^k*prod(1:k)))*diff((x.^2-1)^k,k);% 式 3.65
end

  拟合完成勒让德多项式后，需要求解勒让德多项式的极值点， n 阶勒让德多项式共有 n+1 个极值点。代码如下：

% 求解出 n 阶 Legendre 拟合多项式的极值点
T_diff(x)= diff(T(x),x);  % 式 3.91
x_root_q=solve(T_diff==0);  % 式 3.92

  在获得到极值点后，需要对具有虚部以及超出标准区间范围的标准极值点进行归一化的操作，进行归一化的公式和 matlab 代码如下：

xq,r(i)={1,Im⁡(xq,r(i))≠0−1,Im⁡(xq,r(i))=0 and ∣Re⁡(xq,r(i))∣>1,r=1,2,…,n−1yjk,i,Im⁡(xq,r(i))=0 and ∣Re⁡(xq,r(i))∣≤1x_{q, r}^{(i)}=\left\{\begin{array}{l} 1, \quad \operatorname{Im}\left(x_{q, r}^{(i)}\right) \neq 0 \\ -1, \quad \operatorname{Im}\left(x_{q, r}^{(i)}\right)=0 \text { and }\left|\operatorname{Re}\left(x_{q, r}^{(i)}\right)\right|>1, \quad r=1,2, \ldots, n-1 \\ y_{j}^{k, i}, \operatorname{Im}\left(x_{q, r}^{(i)}\right)=0 \text { and }\left|\operatorname{Re}\left(x_{q, r}^{(i)}\right)\right| \leq 1 \end{array}\right. 

% x_root 为第 q 个区间参数的极值点集合，最终的所有区间参数的极值点由这些集合中的最大和最小值组成
x_root_q(imag(x_root_q)~=0)=1;% 式 3.93 虚数等于 1
x_root_q(abs(x_root_q)>1)=-1;% 式 3.93 超出的区间的数等于 - 1
x_root{q}=[double(x_root_q)',-1,1];% 存储所有极值点

  在完成极值点的标准化后，需要将标准极值点和边界值 -1 和 1 一起使用拟合完成的勒让德多项式进行运算，分别得到对应勒让德多项式最大 / 最小的参考点（某个极值点或 - 1 或 1），在找到响应最大 / 最小参考点后，根据当前参数区间从 [-1, 1] 映射到响应的参数区间，得到对应最大 / 最小目标函数值的真实输入点。响应的 matlab 代码如下：

T_x_root=T(x_root{q}); % 式 3.94 找出找出第 q 个分量的最大最小值点
x_root_min=x_root{q}(find(T_x_root==min(T_x_root)));% 式 3.95 可能会出现多个最小值，此时必为区间边界，取其中任意一个即可
x_root_max=x_root{q}(find(T_x_root==max(T_x_root)));% 式 3.95 可能会出现多个最大值，此时必为区间边界，取其中任意一个即可
p_min(i)=double(vpa(pC_All(i)+pR_All(i).*x_root_min(1)));% 式 3.96+3.97
p_max(i)=double(vpa(pC_All(i)+pR_All(i).*x_root_max(1))); % 式 3.96+3.97

  在对所有输入维度进行逐维运算后，需要将结果进行保存，响应的 matlab 代码如下：

u_min=MyFun(p_min);% 式 3.98 此时，解中包含其他维度的值
u_max=MyFun(p_max);% 式 3.98 此时，解中包含其他维度的值
U(q,:)=[u_min(q),u_max(q)]% 最终要求的 U

  至此 U 中即包含了使用区间逐维得到的区间运算结果。其中第ii 行第11 列代表了第ii 个输出维度的最小值，其中第ii 行第22 列代表了第ii 个输出维度的最大值。

# 三、区间逐维运算的全部代码

XMH_IDWA2.m（区间逐维运算的主函数）
function U=XMH_IDWA2(MyFun,pMin,pMax,m,N,n,s)
pC_All=(pMin+pMax)/2;  % 中值
pR_All=(pMax-pMin)/2;  % 半径
%% 计算样本点对应的函数值
for i=1:m  % 对每个优化参数进行一次区间逐维分析
    pC=pC_All(i);  % 第 i 个参数的中值
    pR=pR_All(i);  % 第 i 个参数的半径
    x_gauss=Gauss(s);  % 式 3.82 获取高斯积分点 有 s+1 个
    p_gauss=pC+pR.*x_gauss;  % 式 3.83 + 式 3.84 获取第 i 个参数的真实采样点
    
    % 式 3.85 + 式 3.86
    S_p_gauss=repmat(pC_All, s+1, 1);
    for j=1:s+1
        S_p_gauss(j,i)=p_gauss(j);
    end
    
    % 式 3.87 + 式 3.88
    for j=1:s+1
        u_all=MyFun(S_p_gauss(j,:));% 此时，解中包含其他维度的值
        cd(j+(s+1)*(i-1),:)=u_all;% 保存所有真实计算的结果 列是响应向量 行是
    end
end
%% 逐维
for q=1:N % 目标函数的维度
    for i=1:m % 对每个优化参数进行一次区间逐维分析
        u=cd((s+1)*(i-1)+1:(s+1)*i,q);% 对应参数和对应目标函数的函数关系
        % 求解 n 阶 Legendre 拟合多项式
        syms x
        for k=0:n
            % 求 A_Tq_k_1
            A_Tq_k_1=0;
            for j=1:s+1
                if k==0
                    L_k_x=1;
                else
                    L_k(x)=Legendre_k(k,x);
                    L_k_x=L_k(x_gauss(j));
                end
                L_s_diff(x)=diff(Legendre_k(s+1,x),1);
                A_Tq_k_1=A_Tq_k_1+((2*k+1)*u(j)*L_k_x)/((1-x_gauss(j)^2)*vpa((L_s_diff(x_gauss(j)))^2));% Lk (x) 的系数 公式 3.90
            end
            % 利用 k 和 A_Tq_k_1 求解 T (x)（每一次都累加，共累加 n 次，n 阶平方逼近）
            if k==0
                T(x)=A_Tq_k_1*Legendre_k(k,x);
            else
                T(x)=T(x)+A_Tq_k_1*Legendre_k(k,x);
            end
        end
        % 求解出 n 阶 Legendre 拟合多项式的极值点
        T_diff(x)= diff(T(x),x);  % 式 3.91
        x_root_q=solve(T_diff==0);  % 式 3.92
        
        % x_root 为第 q 个区间参数的极值点集合，最终的所有区间参数的极值点由这些集合中的最大和最小值组成
        x_root_q(imag(x_root_q)~=0)=1;% 式 3.93 虚数等于 1
        x_root_q(abs(x_root_q)>1)=-1;% 式 3.93 超出的区间的数等于 - 1
        x_root{q}=[double(x_root_q)',-1,1];% 存储所有极值点
        T_x_root=T(x_root{q}); % 式 3.94 找出找出第 q 个分量的最大最小值点
        x_root_min=x_root{q}(find(T_x_root==min(T_x_root)));% 式 3.95 可能会出现多个最小值，此时必为区间边界，取其中任意一个即可
        x_root_max=x_root{q}(find(T_x_root==max(T_x_root)));% 式 3.95 可能会出现多个最大值，此时必为区间边界，取其中任意一个即可
        
        p_min(i)=double(vpa(pC_All(i)+pR_All(i).*x_root_min(1)));% 式 3.96+3.97
        p_max(i)=double(vpa(pC_All(i)+pR_All(i).*x_root_max(1))); % 式 3.96+3.97
        
    end
    u_min=MyFun(p_min);% 式 3.98 此时，解中包含其他维度的值
    u_max=MyFun(p_max);% 式 3.98 此时，解中包含其他维度的值
    
    U(q,:)=[u_min(q),u_max(q)]% 最终要求的 U
    
end
end
Gauss.m（高斯拟合的函数）
function y= Gauss(n)
    % gauss-legendre 求积法点集的求解程序 论文里原代码
    if n==0
        y=0;
    else
        beta = 0.5./sqrt(1-(2*(1:n)).^(-2));% B = sqrt (X) 返回数组 X 的每个元素的平方根
        T = diag(beta,1)+diag(beta,-1);% 取 beta 中的元素为上方第 1 条对角线和下方第一条对角线
        [~,D]=eig(T);% eig 为求特征值和特征向量
        x=diag(D);% 取特征值
        [x,~]=sort(x);
        y=x;
    end
end
Legendre_k.m（用于拟合勒让德多项式的函数）
function L_k=Legendre_k(k,x)
    syms x
    L_k=1./double((2.^k*prod(1:k)))*diff((x.^2-1)^k,k);% 式 3.65
end
KUF.m（算例使用的函数）
function u=KUF(p)
%%%% 计算响应值 %%%%%%
K=[p,1/p;(p-2)^2,1];
F=[2*p;p^2];
u=pinv(K)*F;
%%%%%%%%%%%%%%
end
main.m（算例主函数）
clear all
PL=[1.5];% 区间变量下界
PU=[2.5];% 区间变量上界
PM=length(PL);% 区间变量维数
PN=2;% 目标函数的维度
MyFun=@KUF;% 函数
Ln=3;% 几阶平方逼近
Ls=3;% 高斯积分点的个数
Mn=10000;% 蒙特卡洛评估次数
U1=XMH_IDWA2(MyFun,PL,PU,PM,PN,Ln,Ls); % 区间逐维
disp('区间逐维方法结果如下：');% 显示结果
disp(U1);% 显示结果
U2=Monte_Carlo(MyFun,Mn,PL,PU,PM,PN); % 蒙特卡洛
disp('蒙特卡洛方法结果如下：');% 显示结果
disp(U2);% 显示结果

      

          
          
      

  


   


    

  
    
      
    
    更新于
    
  
  
      
        
      
      阅读次数
      
      次