您所在位置：主题：[讨论]我的这个程序总是不出优化结果，请高手指点
qinfeirong
[专家分：0] 发布于 2010-06-08 17:47:00
有么有人能帮帮我，我的未知数有七个，怎么也出不了结果
%七杆机构优化(主函数) linkopt_main
function linkopt_main
clc;
clear all;
global L1
L1=2400; %L1为l7=l9
x0=[115;865;862;800;988;1624;117]; %优化变量赋初值,对应[l1,l3,l5,l6,l10,l11,θ10]
lb=[112,800,800,700,900,1400,110]
ub=[140,900,900,900,1200,1800,130]
global x;
options=optimset('LargeScale','on','TolFun',1e-12);
[x,feval,exitflag,output]=fmincon(@link_objfun,x0,[],[],[],[],lb,ub,@link_confun,options)
%目标函数
function f=link_objfun(x)
global L1
global AC
global r
global L8
global a
global t
global GL
global tz
global Mm
global by
global v
global b
global w
global m
global f
L1=2400; %L1为l7=l9
AC=x(1)*cos(1/180*(270-x(7))*pi)+(x(1)^2*cos(1/180*(270-x(7))*pi)^2+x(2)^2-x(1)^2)^(1/2) %根据有关AC的方程x(2)^2-x(1)^2-AC^2+2*x(1)*AC*cos((270-x(7))/180*pi)=0反解出
r(1)=atan((AC-x(5))/x(6)) ; %为角度γ1
DH=sqrt((AC-x(5))^2+x(6)^2)
r(2)=acos((DH^2+x(3)^2-x(4)^2)/(2*DH*x(3))) ; %为角度γ2
r(3)=pi-(r(2)-r(1)) %为角度γ3
L8=sqrt(x(3)^2+L1^2-2*x(3)*L1*cos(r(3))) %L8是一个变量，但是它可以有其他变量推出，所以不算在优化变量内
a=pi-x(7)/180*pi %a指初始角度θ20
t=acos(((x(3))^2+L1^2-L8^2)/(2*(x(3))*L1)); %为角CDG
GL=sqrt(x(3)^2+(L1^2)/4-x(3)*L1*cos(t));
tz=asin((L1*sin(t))/(2*GL)); %为角LGD
Mm=[cos(t),sin(t),GL*cos(tz);-sin(t),cos(t),GL*sin(tz);0,0,1]; %为坐标变换矩阵之一
by=acos((x(3)^2+(L8)^2-L1^2)/(2*x(3)*L8)) ; %角CGD
f=0;
g=0;
h=0;
for i=75:107 %将杆1一周角度划分为180等分，循环的初始却是在i=75时。该推导是由仿真决定的
v=117/180*pi+2/180*pi*i %曲柄转角各等分点的值
%求解相关角度
qiujiejiaodu;
%求解角速度
jiaosudu;
%求解动瞬心
jieshunxin;
%进行坐标变换
Mp=[cos(b(5)),-sin(b(5)),x(4)*cos(b(1));sin(b(5)),cos(b(5)),x(4)*sin(b(1));0,0,1]; %为坐标变换矩阵之一,随角度变化而变化
u=Mp*Mm*[m(1);m(2);1];
g=g+(u(2)+897)^2;
h=h+((sqrt(m(1)^2+(m(2)-71079)^2)-72000))^2;
end %for循环结束
f=f+sqrt(g/33)+sqrt(h/33);
%约束函数
function [c ceq]=link_confun(x,L1,L8,a,Mm,b,t,by)
global L1
global x
global L8
global a
global Mm
global b
global t
global by %给f赋初值
c=[205-x(1)-x(2)+x(2)*sqrt(1-(x(1)/x(2))^2*(cos(x(7)/180*pi))^2)-x(1)*sin(x(7)/180*pi)];
ceq=[];
下面是一些调用小函数
function b=qiujiejiaodu
b0=[-3.5;1.1;2;1.05;2.5]; %初始点
options=optimset('Display','off'); %不显示输出信息
global b
b=fsolve(@jd,b0,options)
%这个求解角度的函数，输出的是f值，是一个五行一列的矩阵。
function f=jd(b,v,a,L8,L1,x,by) %求解角度是为了后面求解角速度，而求解角速度是为了在求解速度瞬心的矩阵中加以应用
global v
global L8
global a
global L1
global x
global by
%global b
f = [x(4)*cos(b(1))+x(3)*cos(b(5))+x(2)*cos(b(4))+x(6)-x(1)*cos(b(2));
x(4)*sin(b(1))+x(3)*sin(b(5))+x(2)*sin(b(4))-x(5)-x(1)*sin(b(2));
x(4)*cos(b(1))+L8*cos(b(5)+by)+x(2)*cos(b(3))+L1+x(6)-x(1)*cos(v);
x(4)*sin(b(1))+L8*sin(b(5)+by)+x(2)*sin(b(3))-x(5)-x(1)*sin(v);
v-(x(7)/180*pi-a)-b(2)]; %b(1)代表θ6，其余各b对应相应θ
function w=jiaosudu
w0=[1;1;1;1;1]; %初始点
options=optimset('Display','off'); %不显示输出信息
global w
w = fsolve(@jsd,w0,options)
function f=jsd(w,L1,L8,x,b,v,by)
global by
global b
global v
global L8
global L1
global x
global cw
cw=2; %给l1的旋转赋角速度初值
f=[(-x(4))*sin(b(1))*w(1)-x(3)*sin(b(5))*w(5)-x(2)*sin(b(4))*w(4)+x(1)*sin(b(2))*w(2);
x(4)*cos(b(1))*w(1)+x(3)*cos(b(5))*w(5)+x(2)*cos(b(4))*w(4)-x(1)*cos(b(2))*w(2);
(-x(4))*sin(b(1))*w(1)-L8*sin(b(5)+by)*w(5)-x(2)*sin(b(3))*w(3)+x(1)*sin(v)*cw;
x(4)*cos(b(1))*w(1)+L8*cos(b(5)+by)*w(5)+x(2)*cos(b(3))*w(3)-x(1)*cos(v)*cw;
w(2)-cw]; %b(1)代表θ6，其余各b对应相应θ,w(1)代表ω（6），cw代表ω（1）
function m=jieshunxin %求解动瞬心
m0=[0;0]; %初始点
options=optimset('Display','off'); %不显示输出信息
global m
m=fsolve(@jsx,m0,options)
function f=jsx(m,b,w,x,Mm)
global b
global w
global x
global Mm
Ms=[-sin(b(5))*w(5),-cos(b(5))*w(5),(-x(4))*sin(b(1))*w(1);cos(b(5))*w(5),-sin(b(5))*w(5),(x(4))*cos(b(1))*w(1);0,0,0];
f=[Ms*Mm*[m(1);m(2);1]];
将每一个函数都保存了，运行主函数linkopt_main就能知道错误提示了，奇怪的是为什么总是说x(4)找不到。
有哪位高手能解开这个问题，小女子感激不尽，这是连老师也无奈的程序了
%七杆机构优化(主函数) linkopt_main
function linkopt_main
clc;
clear all;
global L1
L1=2400; %L1为l7=l9
x0=[115;865;862;800;988;1624;117]; %优化变量赋初值,对应[l1,l3,l5,l6,l10,l11,θ10]
lb=[112,800,800,700,900,1400,110]
ub=[140,900,900,900,1200,1800,130]
global x;
options=optimset('LargeScale','on','TolFun',1e-12);
[x,feval,exitflag,output]=fmincon(@link_objfun,x0,[],[],[],[],lb,ub,@link_confun,options)
%目标函数
function f=link_objfun(x)
global L1
global AC
global r
global L8
global a
global t
global GL
global tz
global Mm
global by
global v
global b
global w
global m
global f
L1=2400; %L1为l7=l9
AC=x(1)*cos(1/180*(270-x(7))*pi)+(x(1)^2*cos(1/180*(270-x(7))*pi)^2+x(2)^2-x(1)^2)^(1/2) %根据有关AC的方程x(2)^2-x(1)^2-AC^2+2*x(1)*AC*cos((270-x(7))/180*pi)=0反解出
r(1)=atan((AC-x(5))/x(6)) ; %为角度γ1
DH=sqrt((AC-x(5))^2+x(6)^2)
r(2)=acos((DH^2+x(3)^2-x(4)^2)/(2*DH*x(3))) ; %为角度γ2
r(3)=pi-(r(2)-r(1)) %为角度γ3
L8=sqrt(x(3)^2+L1^2-2*x(3)*L1*cos(r(3))) %L8是一个变量，但是它可以有其他变量推出，所以不算在优化变量内
a=pi-x(7)/180*pi %a指初始角度θ20
t=acos(((x(3))^2+L1^2-L8^2)/(2*(x(3))*L1)); %为角CDG
GL=sqrt(x(3)^2+(L1^2)/4-x(3)*L1*cos(t));
tz=asin((L1*sin(t))/(2*GL)); %为角LGD
Mm=[cos(t),sin(t),GL*cos(tz);-sin(t),cos(t),GL*sin(tz);0,0,1]; %为坐标变换矩阵之一
by=acos((x(3)^2+(L8)^2-L1^2)/(2*x(3)*L8)) ; %角CGD
f=0;
g=0;
h=0;
for i=75:107 %将杆1一周角度划分为180等分，循环的初始却是在i=75时。该推导是由仿真决定的
v=117/180*pi+2/180*pi*i %曲柄转角各等分点的值
%求解相关角度
qiujiejiaodu;
%求解角速度
jiaosudu;
%求解动瞬心
jieshunxin;
%进行坐标变换
Mp=[cos(b(5)),-sin(b(5)),x(4)*cos(b(1));sin(b(5)),cos(b(5)),x(4)*sin(b(1));0,0,1]; %为坐标变换矩阵之一,随角度变化而变化
u=Mp*Mm*[m(1);m(2);1];
g=g+(u(2)+897)^2;
h=h+((sqrt(m(1)^2+(m(2)-71079)^2)-72000))^2;
end %for循环结束
f=f+sqrt(g/33)+sqrt(h/33);
%约束函数
function [c ceq]=link_confun(x,L1,L8,a,Mm,b,t,by)
global L1
global x
global L8
global a
global Mm
global b
global t
global by %给f赋初值
c=[205-x(1)-x(2)+x(2)*sqrt(1-(x(1)/x(2))^2*(cos(x(7)/180*pi))^2)-x(1)*sin(x(7)/180*pi)];
ceq=[];
下面是一些调用小函数
function b=qiujiejiaodu
b0=[-3.5;1.1;2;1.05;2.5]; %初始点
options=optimset('Display','off'); %不显示输出信息
global b
b=fsolve(@jd,b0,options)
%这个求解角度的函数，输出的是f值，是一个五行一列的矩阵。
function f=jd(b,v,a,L8,L1,x,by) %求解角度是为了后面求解角速度，而求解角速度是为了在求解速度瞬心的矩阵中加以应用
global v
global L8
global a
global L1
global x
global by
%global b
f = [x(4)*cos(b(1))+x(3)*cos(b(5))+x(2)*cos(b(4))+x(6)-x(1)*cos(b(2));
x(4)*sin(b(1))+x(3)*sin(b(5))+x(2)*sin(b(4))-x(5)-x(1)*sin(b(2));
x(4)*cos(b(1))+L8*cos(b(5)+by)+x(2)*cos(b(3))+L1+x(6)-x(1)*cos(v);
x(4)*sin(b(1))+L8*sin(b(5)+by)+x(2)*sin(b(3))-x(5)-x(1)*sin(v);
v-(x(7)/180*pi-a)-b(2)]; %b(1)代表θ6，其余各b对应相应θ
function w=jiaosudu
w0=[1;1;1;1;1]; %初始点
options=optimset('Display','off'); %不显示输出信息
global w
w = fsolve(@jsd,w0,options)
function f=jsd(w,L1,L8,x,b,v,by)
global by
global b
global v
global L8
global L1
global x
global cw
cw=2; %给l1的旋转赋角速度初值
f=[(-x(4))*sin(b(1))*w(1)-x(3)*sin(b(5))*w(5)-x(2)*sin(b(4))*w(4)+x(1)*sin(b(2))*w(2);
x(4)*cos(b(1))*w(1)+x(3)*cos(b(5))*w(5)+x(2)*cos(b(4))*w(4)-x(1)*cos(b(2))*w(2);
(-x(4))*sin(b(1))*w(1)-L8*sin(b(5)+by)*w(5)-x(2)*sin(b(3))*w(3)+x(1)*sin(v)*cw;
x(4)*cos(b(1))*w(1)+L8*cos(b(5)+by)*w(5)+x(2)*cos(b(3))*w(3)-x(1)*cos(v)*cw;
w(2)-cw]; %b(1)代表θ6，其余各b对应相应θ,w(1)代表ω（6），cw代表ω（1）
function m=jieshunxin %求解动瞬心
m0=[0;0]; %初始点
options=optimset('Display','off'); %不显示输出信息
global m
m=fsolve(@jsx,m0,options)
function f=jsx(m,b,w,x,Mm)
global b
global w
global x
global Mm
Ms=[-sin(b(5))*w(5),-cos(b(5))*w(5),(-x(4))*sin(b(1))*w(1);cos(b(5))*w(5),-sin(b(5))*w(5),(x(4))*cos(b(1))*w(1);0,0,0];
f=[Ms*Mm*[m(1);m(2);1]];
将每一个函数都保存了，运行主函数linkopt_main就能知道错误提示了，奇怪的是为什么总是说x(4)找不到。
有哪位高手能解开这个问题，小女子感激不尽，这是连老师也无奈的程序了