随机方向法matlab源代码

问题以机械优化设计（哈工大第五版课本例6-2）为例， 求解如下图所示的约束优化问题。

源代码包含三个matlab的m文件，分别为目标函数、约束函数及随机方向法主文件。

（1）目标函数fun.m

function f=fun(x)
f=x(1)^2+x(2);

(2)约束函数gcon.m

function g=gcon(x)
g=[
    x(1)^2+x(2)^2-9
    x(1)+x(2)-1
];

(3)随机方向法主文件RandomDirectionMethod.m

clc,clear
tic  %计时开始
%% 初始量赋值（目标函数、约束条件单独写到m文件中）
n=2;    %自变量个数
a0=0.1; %试验步长
epsilon=1e-6;   %收敛精度
k=8;    %随机搜索方向个数
a=-3;   %自变量下限值
b=3;    %自变量上限值
%%

%% 选择初始点
while 1
    x0=a+(b-a)*rand(n,1);  %产生一个随机初始点，由约束1得自变量取值区间为（-3，3）
    g=gcon(x0);   %计算产生的随机点约束函数值
    if all(g<=0)
        break;
    end
end
f0=fun(x0);  %计算初始随机可行点函数值
% x0,f0  %调试是否生成初始可行点x0
%%

while 1
    %% 产生可行搜索方向
    a0=0.1;
    r=zeros(n,k);   %存放k个随机方向
    e=zeros(n,k);   %存放k个单位化随机方向
    while 1
        j=1;    %用来记录沿随机方向移动后可行点数目的计数变量
        for i=1:k
          r(:,i)=rands(n,1);    %产生一个坐标（-1，1）之间的随机方向
          e(:,i)=r(:,i)/norm(r(:,i));   %单位化随机方向
          xs=x0+a0*e(:,i);  %计算沿随机方向移动后的自变量点
          gz=gcon(xs);          
          if all(gz<=0)   
            fs=fun(xs); %计算移动后可行点的目标函数值
            x(:,j)=xs;  %记录沿随机方向移动后的可行点
            fz(:,j)=fs; %记录沿随机方向移动后的可行点目标函数值
            j=j+1;             
          else
             continue;
          end
        end
        xf=[x;fz];  %xf为三行矩阵，前两行代表移动后可行点，第三行代表移动后可行点的目标函数值
        [B,ind]=sort(xf(n+1,:));    %升序排列移动后可行点的目标函数值
        fl=B(1);    %获取最小的移动后可行点的目标函数值
        if fl<f0
            break;
        else
            a0=0.9*a0;  %若的所有可行点目标函数值都大于初始值f0,则缩小步长，重新生成随机方向试探
        end
    end
    xl=x(:,ind(1)); %找到沿随机方向移动后目标函数值最小且小于f0的可行点
    d=xl-x0;    %找到一个可行搜索方向
    % xl,fl,f0,d  %调试是否找到一个可行方向
    f0=fl;  %将初始点移动到试探可行点
    x0=xl;
    %%
    
    %% 沿可行方向采用加速步长进行搜索
    while 1
        a0=1.3*a0;  %步长加速
        xl=x0+a0*d; %计算沿可行方向移动到的下一点
        fl=fun(xl); %计算移动后的目标函数值
        g=gcon(xl); %计算移动后的约束函数值
        if(all(g<=0) && fl<f0)
            x0=xl;  %若移动的点满足约束条件且函数值下降，则将x0点移动到xl
            f0=fl;
        else
            break;
        end        
    end
    xn=x0-a0*d/1.3; %计算满足约束条件且函数值下降的极限点x0的前一点xn
    fl=fun(xn);     %计算满足约束条件且函数值下降的极限点x0的前一点xn的函数值，以计算收敛条件
    % xn,fl
    % x0,f0 %调试沿可行方向加速搜索后的可行点
    %%
    
    if abs((f0-fl)/f0)<epsilon  %判断是否达到收敛精度
        break;
    end
end
xe=x0;
fe=f0;
xe,fe  %输出最优点坐标及函数值
toc  %计时结束

程序运行结果如下：

 由于初始点是程序随机生成的，故每次运行结果都有细微差别。可以看到运行结果和理论解（0，-3），3十分接近，通过设置改小收敛精度epsilon的值，可以提高计算精度。

为防丢失，源程序文件夹放在了下面百度网盘链接：

链接：https://pan.baidu.com/s/1KQ8uWSGcZd7gYCBdX0xTCw 
提取码：gsse