本帖最后由 wwwjjj7008 于 2018-3-21 10:24 编辑
问题描述:从四个参数的给定范围r1[200 300]、r2[50 100]、l1[150 200]、l2[150 250]中,任取一组数,代入:
h = sqrt(l2*l2-(3*r1/2-3*l1*cos(t)/2-90.9).^2 -(3^0.5*r1/2-3^0.5*l1*cos(t)/2-119.1).^2);
x = 0.39*r1-0.16*r2-0.39*l1*cos(t)-0.54*h+8.24 ;
z = 1.39*r1 - 1.16*r2 - 1.39 * l1 * cos(t) + l1*sin(t) +0.46*h+14.41;
其中t的变化范围是:50*pi/180~57*pi/180
得到xz曲线,求曲线关于坐标轴所围成面积最大时的那组参数,并输出面积和最大时的图像。
下面是我用fmincon,globalsearch和multistart分别的结果:
GetS.m:
function S = GetS(data)
% 待优化函数
% 输出S 为待优化的面积 z(x)与x轴围成的面积
r1 = data(1);r2 = data(2);l1 = data(3);l2 = data(4);
syms t;
h = sqrt(l2*l2-(3*r1/2-3*l1*cos(t)/2-90.9).^2 -(3^0.5*r1/2-3^0.5*l1*cos(t)/2-119.1).^2);
x = 0.39*r1-0.16*r2-0.39*l1*cos(t)-0.54*h+8.24 ;
z = 1.39*r1 - 1.16*r2 - 1.39 * l1 * cos(t) + l1*sin(t) +0.46*h+14.41;
% x(t)的一阶导数
x_ = diff(x);
del = 100;
t = linspace(50*pi/180,57*pi/180,del);
% 将t代入上述表达式中 求解对应的x,z
xda = subs(x);
zda = subs(z);
x_da = subs(x_);
digits(10);
% 将数据格式转化为double
zda = double(vpa(zda));
xda = double(vpa(xda));
x_da = double(vpa(x_da));
% 忽略虚数解
zda = zda(abs(imag(zda))<=1e-6);
xda = xda(abs(imag(xda))<=1e-6);
x_da = x_da(abs(imag(x_da))<=1e-6);
% 参数方程的积分方法
deltat = pi/(2*del);
% 积分 原理有附图解释
S = sum(x_da.*zda*deltat);
% 用trapz也可以
% S = trapz(xda,zda);
% 为了优化 取负
S = -S;
end
plotxz.m:
function y = plotxz( data )
% 输入函数data为数组 分别对应 r1 r2 l1 l2
% 函数功能:画函数图像z(x)
% data分别对应r1 r2 l1 l2
r1 = data(1);r2 = data(2);l1 = data(3);l2 = data(4);
% 定义符号变量 t
syms t;
h = sqrt(l2*l2-(3*r1/2-3*l1*cos(t)/2-90.9).^2 -(3^0.5*r1/2-3^0.5*l1*cos(t)/2-119.1).^2);
x = 0.39*r1-0.16*r2-0.39*l1*cos(t)-0.54*h+8.24 ;
z = 1.39*r1 - 1.16*r2 - 1.39 * l1 * cos(t) + l1*sin(t) +0.46*h+14.41;
% t 生成100个50*pi/180~57*pi/180之间的数
t = linspace(50*pi/180,57*pi/180,100);
% 将t代入上述表达式中 求解对应的x,z
xda = subs(x);
zda = subs(z);
% 将数据格式转化为double
zda = double(vpa(zda));
xda = double(vpa(xda));
% 忽略虚数解
zda = zda(abs(imag(zda))<=1e-6);
xda = xda(abs(imag(xda))<=1e-6);
figure
plot(xda,zda)
grid on
y=1;
end
%下面是主函数
% 设置优化约束条件
options.TolCon = 1e-15;
options.TolX = 1e-15;
options.TolFun = 1e-15;
% 设置迭代优化代数:500
options = gaoptimset(options, 'Generations', 500);
%上边界
LB = [200,50,150,150];
% 上边界
UB = [300,100,200,250];
% 优化初始值
x0 = [200,50,180,180];
关键是下面用三种不同的优化主函数,得到的结果不同,而且对初始值非常敏感:
(1)用fmincon优化的主函数
main.m:
options = optimset('LargeScale','off');
% 优化函数 fmincon 待优化函数为 GetS函数
[ANSW,FVAL] = fmincon(@(data)GetS(data),x0,[],[],[],[],LB,UB,[],options);
% 优化结果r1 r2 l1 l2 存放在ANSW数组里面
plotxz( ANSW );
% 优化求解的最大面积
S=-FVAL
结果1.png (13.45 KB, 下载次数: 1)
2018-3-21 10:22 上传
结果2.png (13.25 KB, 下载次数: 1)
2018-3-21 10:22 上传
(2)用globalsearch优化的主函数
main.m:
f = @(data)GetS(data);
problem = createOptimProblem('fmincon','objective',f,'x0',x0,'lb',LB,'ub',UB,'options',optimset('Algorithm','SQP','DISP','none'));
gs = GlobalSearch;
xgs = run(gs,problem);
% 优化结果r1 r2 l1 l2 存放在xgs数组里面
plotxz( xgs );
xgs
S = -GetS(xgs)
结果1.png (7.58 KB, 下载次数: 1)
2018-3-21 10:23 上传
结果2.png (12.35 KB, 下载次数: 1)
2018-3-21 10:23 上传
(3)用multistart优化的主函数
main.m:
f = @(data)GetS(data);
problem = createOptimProblem('fmincon','objective',f,'x0',x0,'lb',LB,'ub',UB,'options',optimset('Algorithm','SQP','DISP','none'));
ms = MultiStart;
tic
xgs_m = run(ms,problem,50);
toc
% 优化结果r1 r2 l1 l2 存放在xgs_m数组里面
plotxz( xgs_m );
xgs_m
S = -GetS(xgs_m)
结果1.png (14.02 KB, 下载次数: 0)
2018-3-21 10:23 上传
结果2.png (11.04 KB, 下载次数: 2)
2018-3-21 09:55 上传
求大神指点迷津,我这到底哪里出错了?理论上globalsearch和multistart的结果应该会比fmincon的结果更优不是吗?