【求助】悬置MATLAB优化

LAN博基尼

大家好，我现在有一振动模型，具体算法如下：

1. clear
   
2. clc
   
3. %动力总成质量
   
4. m=323.98;%kg
   
5. %惯性积
   
6. Jx=25.706;Jy=13.607;Jz=23.014;Jxy=-3.029;Jyz=3.359;Jzx=-0.876;%kg*m*m
   
7. %各个悬置点坐标
   
8. Xz=72.44;Yz=-431.92;Zz=148.82;
   
9. Xy=4.34;Yy=452.98;Zy=250.22;
   
10. Xh=199.64;Yh=-7.72;Zh=-245.93;
    
11. %mm在下面将位移坐标变成m
    
12. %各个悬置的三向刚度
    
13. Kzx=380e3*1.5;Kzy=72e3*1.5;Kzz=220e3*1.5;
    
14. Kyx=380e3*1.5;Kyy=72e3*1.5;Kyz=220e3*1.5;
    
15. Khx=213e3*1.5;Khy=10e3*1.5;Khz=10e3*1.5;%N/m
    
16. %各个悬置的安装角度
    
17. %左悬置安装角度
    
18. ALFAzu=0;BETAzu=90*pi/180;GAMAzu=90*pi/180;
    
19. ALFAzv=90*pi/180;BETAzv=0;GAMAzv=90*pi/180;
    
20. ALFAzw=90*pi/180;BETAzw=90*pi/180;GAMAzw=0;
    
21. %右悬置安装角度
    
22. ALFAyu=0;BETAyu=90*pi/180;GAMAyu=90*pi/180;
    
23. ALFAyv=90*pi/180;BETAyv=0;GAMAyv=90*pi/180;
    
24. ALFAyw=90*pi/180;BETAyw=90*pi/180;GAMAyw=0;
    
25. %后悬置安装角度
    
26. ALFAhu=8*pi/180;BETAhu=90*pi/180;GAMAhu=82*pi/180;
    
27. ALFAhv=90*pi/180;BETAhv=0;GAMAhv=90*pi/180;
    
28. ALFAhw=98*pi/180;BETAhw=90*pi/180;GAMAhw=8*pi/180;
    
29. %安装角矩阵Tj
    
30. Th=[cos(ALFAhu),cos(BETAhu),cos(GAMAhu);cos(ALFAhv),cos(BETAhv),cos(GAMAhv);cos(ALFAhw),cos(BETAhw),cos(GAMAhw)];
    
31. Tz=[cos(ALFAzu),cos(BETAzu),cos(GAMAzu);cos(ALFAzv),cos(BETAzv),cos(GAMAzv);cos(ALFAzw),cos(BETAzw),cos(GAMAzw)];
    
32. Ty=[cos(ALFAyu),cos(BETAyu),cos(GAMAyu);cos(ALFAyv),cos(BETAyv),cos(GAMAyv);cos(ALFAyw),cos(BETAyw),cos(GAMAyw)];
    
33. %参数输入完毕
    
34. %%%%%%%%%%%====================%%%%%%%%%%
    
35. %%%%%%%%%%%---能量表达式中的Fj和kj---%%%%%%%%%%
    
36. Fh=[1e3,0,0,0,Zh,-Yh;0,1e3,0,-Zh,0,Xh;0,0,1e3,Yh,-Xh,0]*0.001;
    
37. Fz=[1e3,0,0,0,Zz,-Yz;0,1e3,0,-Zz,0,Xz;0,0,1e3,Yz,-Xz,0]*0.001;
    
38. Fy=[1e3,0,0,0,Zy,-Yy;0,1e3,0,-Zy,0,Xy;0,0,1e3,Yy,-Xy,0]*0.001;
    
39. %%%此处将位移坐标转变成m
    
40. kh=[Khx,0,0;0,Khy,0;0,0,Khz];
    
41. kz=[Kzx,0,0;0,Kzy,0;0,0,Kzz];
    
42. ky=[Kyx,0,0;0,Kyy,0;0,0,Kyz];
    
43. %%%%%%%%%%%---振动方程中的M,D和K---%%%%%%%%%%
    
44. M=[m,0,0,0,0,0;0,m,0,0,0,0;0,0,m,0,0,0;0,0,0,Jx,-Jxy,-Jzx;0,0,0,-Jxy,Jy,-Jyz;0,0,0,-Jzx,-Jyz,Jz];
    
45. % D=[1200,0,0,0,0,0;0,1600,0,0,0,0;0,0,2000,0,0,0;0,0,0,60,0,0;0,0,0,0,80,0;0,0,0,0,0,100];
    
46. K=Fh'*Th'*kh*Th*Fh+Fz'*Tz'*kz*Tz*Fz+Fy'*Ty'*ky*Ty*Fy;
    
47. %%%%%%%%%%%%%---下面计算系统固有特性---%%%%%%%%%%%%%
    
48. [v,d]=eig(K,M);
    
49. f=diag(sqrt(d)/2/pi)'
    
50. for j=1:6
    
51. for k=1:6
    
52. for l=1:6
    
53. ENERGY(k,l)=M(k,l)*v(k,j)*v(l,j);
    
54. end
    
55. end
    
56. qq=sum(ENERGY);
    
57. qqt=sum(qq);
    
58. dig=[qq/qqt];
    
59. EGEN(:,j)=dig;
    
60. ep1=max(dig);
    
61. ep2=sum(dig);
    
62. eper(j)=ep1/ep2*100;
    
63. end
    
64. eper
    
65. EGEN=EGEN*100

复制代码

程序前面是参数的输入和一些坐标变换，最后求出来的是系统的固有频率f，能量矩阵EGEN（对应6乘6的矩阵），eper对应EGEN每列最大值。
我现在想提高eper的值，理论上可以通过修改刚度K得到，现在以
object=-sum（max（EGEN））为目标函数：
约束条件为：
（1）程序里面的Kzx、Kzy、Kzz；Kyx、Kyy、Kyz；Khx、Khy、Khz，想把这些刚度值设一个上限和下下限。
（2）还要把求出来的f（i）作为约束条件，设定6个f（i）的范围。
我参考了fmincon、fgoalattain等多个函数，也做过了一些例子，但是对于自己这个问题还是无从下手，希望大家能提供一点思路和帮助，现在急需最后优化的刚度值。
在线等候各位的提示，谢谢，不甚感激。。。
我的邮箱：zs.jordan@163.com

zhouyang664

这个应该到振动或是物理那些论坛去看一下...

bozheng_1984

我也有此问题 ,请联系wliang_auto@hotmail.com.3q

happy

这个问题用fmincon应该可以实现，大概考虑了一下，基本格式应该可以通过下面的格式实现
x = fmincon(@myfun,x0,A,b,Aeq,beq,lb,ub,@mycon)

object=-sum（max（EGEN））
对应函数中的@myfun
约束条件1：对应A,b
约束条件2：对应@mycon

估计你的主要问题应该出现在约束条件2上
具体你可以参考帮助文件中的nonlcon项的相关说明

Happy99

LZ的另一相同意议题的帖
http://forum.vibunion.com/thread-97216-1-1.html

LAN博基尼

嘿嘿，经过几天的调试，已经做出来了。
开始就是约束函数没注意，频率范围设置的时候要包含原始值在内，不然计算不收敛。
谢谢各位。

苍鹰逐梦

LAN博基尼 发表于 2010-11-22 12:35
嘿嘿，经过几天的调试，已经做出来了。
开始就是约束函数没注意，频率范围设置的时候要包含原始值在内，不 ...

您好，我也用matlab中fmincon编了解耦程序，跟你说的差不多，频率设置时候需要注意。但是这样一来，
如果想要优化到想要的频率，那就无法计算了，用这个fmincon就只能增大解耦率。请问你的程序是吗？
谢谢