基于动态径向基函数（DRBF）代理模型的优化策略

基于动态径向基函数（DRBF）代理模型的优化策略

在工程计算中，我们经常遇到需要求解优化问题，尤其在现在十分有前景的机器学习领域，如何快速、高效的求解优化问题，成为机器学习算法是否高效、准确的必要条件。

运用代理模型求解优化问题是优化领域的一个重要思想，该思想是将一个复杂的函数通过采集特征点构造一个相对简单的新的函数，我们称之为代理模型（通常分为插值型与拟合型），再对代理模型进行优化求得最优解。由于该代理的模型只能一定程度上代表真实模型，它到最优解并不能直接代表原目标函数的最优解，所以需要对代理模型进行更新，通常是以一定的判据来重新采点，重新构造代理模型，直到满足收敛条件，便停止对代理模型的更新，所得的最优解即视为原函数的最优解。

本文对采点、代理模型以及代理模型更新策略仅做了简要的介绍，详细内容参见基于计算试验设计与代理模型的飞行器近似优化策略探讨

在众多代理模型中，笔者首先尝试复现的是一种相对简单的代理模型——径向基函数，题中动态一词体现在更新代理模型时，新的采样区间是不断变化的，具体的做法是，以本次求解的最优点为中心，将上一次最优解区间的尺寸同比例缩小为原来的1/ns，在这样的新矩形内重新采集样本点，再以此样本点与原样本点一起构造新的代理模型。

具体的做法参见基于动态径向基函数代理模型的优化策略，在该篇论文中，更新采样区间时所选的中心点为上一次的最优解处，但这样的话存在第一次循环无法调用第零次（因为没有初始化）的问题，在此提供两种思路，一种是进行初始化，如采取初始寻优区间的中心为初始最优解位置；另一种是改为以当前循环的最优解为中心构造采样区间。笔者选择了第二种方法，因为初始最优解的选取在没有先验知识的情况下随意选取似乎有失严谨。同时，龙腾老师也建议笔者采用后一种方法。感谢龙腾老师与李学长、屁屁给予笔者复现时的帮助。

程序是以SC函数（Six-hump camelback function）为目标函数进行代理模型求解的。在运行时有时会遇到矩阵奇异不能计算、并未收敛到最优解便结束循环的情况。出现该情况，一是该算法本身具有缺陷，采样区间仅缩小不扩大会导致4-5次循环后区间缩至很小，这样在构造径向基函数代理模型时要进行矩阵求逆运算时便会遇到矩阵奇异，从而不能求得最优解。这种情况，在当前算法下若想优化，应该适当减小每次迭代中区间的压缩量（即减小论文中的ns值），而对于未寻找到最优解便跳出循环，则是因为收敛条件还不够严格，课修改代码中Is_convergence函数的epsilon值。

下面是matlab代码实现，该SC函数的最优解值应该为-1.0316。

main函数

%% clear data
clc;
clear;
%% initial settings
ga_opts = gaoptimset('Display','off','PopulationSize',30);
nv = 2;
ns = (nv + 1)*(nv + 2)/2;
ns = 5;
lb = [-2 -2];
ub = [2 2];
lbc = lb;
ubc = ub;
c = 4;
count = 0;
bound = ones(2,2);
%% run the first time
X = chooseX(lb,ub,ns,nv);
Y = SC_function(X);
%% forloop
for k = 1:9991
    newfunction = @(x)(invA(X,c)*Y)'*phy2(x,X,c);%construct newfunction
    [xval,fval,exitflag,output]=ga(newfunction,nv,[],[],[],[],lb,ub,[],ga_opts);
    count = count + output.funccount;
    Y_star(k) = SC_function(xval);
    X_star(k,:) = xval;
    if k == 1
        newbound = shrinkspace(X_star(k,:),lbc,ubc,ns);
        lbc = newbound(1,:);
        ubc = newbound(2,:);
        bound = newbound;
        newX = chooseX(newbound(1,:),newbound(2,:),ns,nv);
        X = [X;newX];
        Y = [Y;SC_function(newX)];
        continue;
    end
    if Is_convergence(Y_star(k),Y_star(k-1)) == 1;
        break;
    end
    newbound = shrinkspace(X_star(k,:),lbc,ubc,ns);
    lbc = newbound(1,:);
    ubc = newbound(2,:);
    bound = [bound;newbound];
    newX = chooseX(newbound(1,:),newbound(2,:),ns,nv);
    X = [X;newX];
    Y = [Y;SC_function(newX)];
end 

chooseX

function [ S ] = chooseX( LB,UB,N,D )

S = lhsdesign(N,D,'criterion','maximin');
S = S.*repmat(UB - LB,N,1) + repmat(LB,N,1);
end

invA

function [ output_args ] = invA(x,c)

n = length(x);
A = ones(n,n);
for i = 1:n
    for j = 1:n
        A(i,j) = phy(x(i,:)',x(j,:)',c);
    end
end
output_args = inv(A);
end

Is_convergence

function [ output_args ] = Is_convergence( rk,rk_1 )

epsilon = 0.01;
if (abs((rk - rk_1)/rk_1)) <= epsilon
    output_args = 1;
else
    output_args = 0;
end
end

phy

function [ output_args ] = phy( x1,x2,c )

output_args = ((x1 - x2)'*(x1 - x2)+c^2)^(-0.5);
end

phy2

function [ output_args ] = phy2( x,X,c)

n = length(X);
for i = 1:n
    output_args(i) = phy(x',X(i,:)',c); 
end
output_args = output_args';
end

SC_function

function [ output_args ] = SC_function( X )

output_args = 4*X(:,1).^2 - 2.1 * X(:,1).^4 + 1/3*X(:,1).^6 + ...
    X(:,1).*X(:,2) - 4*X(:,2).^2 + 4*X(:,2).^4;

end

shrinkspace

function [ output_args ] = shrinkspace(x_k,lb_k_1,ub_k_1,ns)

l_k_1 = ub_k_1 - lb_k_1;
lbc = x_k - 1/ns*(l_k_1);
ubc = x_k + 1/ns*(l_k_1);
lbc = max(lbc,[-2 -2]);
ubc = min(ubc,[2 2]);
a = ubc - lbc;
temp = 0.09;
if a(1)<4*temp
    ubc(1) = x_k(1) + 2*temp;
    lbc(1) = x_k(1) - 2*temp;
end
if a(2)<4*temp
    ubc(2) = x_k(2) + 2*temp;
    lbc(2) = x_k(2) - 2*temp;
end
output_args = [lbc;ubc];
end

欢迎大家批评指正。