|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
识别k_youhua使得函数f的值最小
本程序要对参数k_youhua 进行优化识别,他是一个包含十二个值的参数,前七个值要求其全部大于0,后五个值要求大于0小于0.14并且依此递增。
写成数学表达式就是 k_youhua (1-7)>0
0<k_youhua(8)<k_youhua(9)<k_youhua(10)<k_youhua(11)<k_youhua(12)<0.14
最终的结果是要保证函数 f 的值最小
函数迭代的初值为
[[1e10,1e10,1e10,1e10,1e10,1e10,1e10,0.0233,0.0467,0.07 ,0.0933, 0.1167 ]
%M文件如下
function f=CAN_SHU_SHI_BIECAN_shiyan(k_youhua)
L_1 = 0.14; % 结合面部位的长度
L_2 = 0.26 ; % 悬臂梁的长度
L = L_1 + L_2; % 总体的长度
n_1 = 6; % 结合面部分单元数目
n_2 = 7 ; % 悬臂梁部分单元数目
n = n_1 + n_2 ; % 总体的单元数目
x_1 =[0,k_youhua(1,(n_1+2):2*n_1),0.14]; % 结合面部分的x坐标
%对于结合面部分的x坐标应设定一定的约束,如之和为0.14,大于0,小于0.14
x_2 = 0 : L_2 / n_2 : L_2; % 悬臂梁部分的x坐标
x=[x_1 x_2+0.14]; % 总体的x坐标
x(n_1+1)=[];
y = zeros(1,n+1); % 结点的y坐标
% 节点坐标
gNode = [x' y'] ;
% 单元定义
gElement = zeros( n, 3 ) ;
for i=1:n
if i<=n_1
gElement( i, : ) = [ i, i+1, 1 ] ;
else
gElement( i, : ) = [ i, i+1, 2 ] ;
end
end
% 材料性质
h_1 =0.04;
b_1 =0.05;
h = 0.03;
b = 0.05;
I_Y_1 = h_1^3*b_1/12;
I_Y_2 = h^3*b/12;
% 弹性模量 抗弯惯性矩 截面积 密度
gMaterial = [78.5e9, I_Y_1, 0.04*0.05, 7.8e3 % 材料 1
78.5e9, I_Y_2, 0.03*0.05, 7.8e3 ]; % 材料 2
%定义单元体刚度和质量矩阵
syms k me
for i=1:1:4
for j=1:1:4
k(i,j)=0;
me(i,j)=0;
end
end
%定义总体刚度矩阵的大小以便进行组装
[node_number,dummy] = size( gNode ) ;
syms gK gM
for i=1:1:node_number * 2
for j=1:1:node_number * 2
gK(j,i)=0;
gM(j,i)=0;
end
end
%开始设置循环体计算每一个单元的刚度矩阵和质量矩阵并对其进行组装
for ie=1:n
%对于刚度矩阵的建立
E = gMaterial( gElement(ie, 3), 1 ) ;
I = gMaterial( gElement(ie, 3), 2 ) ;
A = gMaterial( gElement(ie, 3), 3 ) ;
xi = gNode( gElement( ie, 1 ), 1 ) ;
yi = gNode( gElement( ie, 1 ), 2 ) ;
xj = gNode( gElement( ie, 2 ), 1 ) ;
yj = gNode( gElement( ie, 2 ), 2 ) ;
L = sqrt( (xj-xi)^2 + (yj-yi)^2 ) ;
k = [ 12*E*I/L^3 6*E*I/L^2 -12*E*I/L^3 6*E*I/L^2
6*E*I/L^2 4*E*I/L -6*E*I/L^2 2*E*I/L
-12*E*I/L^3 -6*E*I/L^2 12*E*I/L^3 -6*E*I/L^2
6*E*I/L^2 2*E*I/L -6*E*I/L^2 4*E*I/L] ;
%对于质量参数的定义
ro = gMaterial( gElement(ie, 3 ), 4 ) ;
me = ro*A*L/420*[156 22*L 54 -13*L
22*L 4*L^2 13*L -3*L^2
54 13*L 156 -22*L
-13*L -3*L^2 -22*L 4*L^2 ] ;
%由上一步得到可以组装的总体刚度矩阵从而对其进行组装
% 把单元刚度和质量矩阵集成到整体刚度矩阵
% ie --- 单元号
for i=1:1:2
for j=1:1:2
for p=1:1:2
for q =1:1:2
m = (i-1)*2+p ;
n = (j-1)*2+q ;
M = (gElement(ie,i)-1)*2+p ;
N = (gElement(ie,j)-1)*2+q ;
gK(M,N) = gK(M,N) + k(m,n) ;
gM(M,N) = gM(M,N) + me(m,n);
end
end
end
end
end
% 为了使刚度矩阵和质量矩阵对称(在计算时可能引入舍入误差)
for i=1:node_number*2
for j=i:1:node_number*2
gK(j,i) = gK(i,j) ;
gM(j,i) = gM(i,j) ;
end
end
%结束刚度矩阵的建立
%加入刚度参数如,即识别刚度参数
gK=eval(gK);
gM=eval(gM);
for i=1:(n_1+1)
gK(i*2-1,i*2-1) = gK(i*2-1,i*2-1) + k_youhua(i);
end
km=inv(gM)*gK;
[P,W2]=eig(km); %求解特征值 和特征向量
f_w=sqrt(W2)/2/pi; %固有频率
f_w=real(f_w);
[heng,lie]=size(f_w);
for i=1:1:heng
pinlv_1(i)=f_w(i,i);
end
pinlv_1=reshape(pinlv_1,1,[]);
pinlv_1=sort(pinlv_1);
%%%%%%%%%%%%%%%%%%%%法向 %%%%%%%%%%%%%%%%%%%%
%保证以下f的值最小
fprintf( '----------- 法向 ----------------\n') ;
f=abs(pinlv_1(1)-225.4)/225.4+abs(pinlv_1(2)-1342)/1342
fprintf( ' 识别频率(法向) \n' ) ;
pinlv_1(1:2)
fprintf( ' 识别刚度(法向) \n' ) ;
k_youhua
fprintf( '--------------------------------------\n\n\n') ;
%M文件到此结束
我用了MATLAB自带的优化函数fmincon进行优化,这个函数什么都用了可是就是不出结果麻烦哪位大侠帮我识别以下k_youhua,用别的方法也行,最用的结果如果 k_youhua(1-7)之和在1e8到1e10之间,可认为结果基本正确,谢谢了
我的Email:wxguang1986@163.com
我急用 谢谢各位了 非常感谢 |
|