声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1631|回复: 3

[综合讨论] fmincon进行优化得不到结果

[复制链接]
发表于 2009-10-14 10:46 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

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
我急用 谢谢各位了 非常感谢
回复
分享到:

使用道具 举报

发表于 2009-10-15 07:23 | 显示全部楼层
和初始值有关 另外可以设置迭代次数等
如果这个函数不可以用 试试优化算法 比如遗传
 楼主| 发表于 2009-10-15 11:07 | 显示全部楼层
迭代次数和初值应该没什么问题 我的设了 不过可以试试遗传算法
发表于 2009-10-23 15:48 | 显示全部楼层
从你给出的初值,可能是没有将数量级差异较大的参数转化为同一数量级。
用fmincon做优化,应该将这些优化求解的参数的数量级用新变量替换为一致,在进行迭代的函数中还原计算。
(例如将优化参数在主程序中除以1e12,在迭代程序中再乘以1e12。)
否则这个优化函数将得不到正确的优化解,你只会看到一些变量有变化,数量级较大的变量几乎不变。

评分

1

查看全部评分

您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-29 22:53 , Processed in 0.063363 second(s), 19 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表