I want to do the optimization but I dont know how to do it ? please help me

Question

Doha Ali 2023-6-27

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1988753-i-want-to-do-the-optimization-but-i-dont-know-how-to-do-it-please-help-me

关闭： Matt J 2023-6-27

clear
clc
%A = eye(4);
A=eye(5);
%b = [80; 320; 50;50];
%UB=[80; 320; 50;50];
b = [80*2; 290*2; 200;200;200];
UB=[80*2; 290*2; 200;200;200];
b = [65*2/2/2; 570/2/2; 200;200;220];
UB=[65*2/2/2; 570/2/2; 200;200;220];
%b = [107; 453; 94;163;200];
%UB=[107; 453; 94;163;200];
%lb = [50;200;30;30];
%LB=[50;200;30;30];
lb = [50;220;50;80;50];
LB=[50;220;50;80;50];
lb = [48*2/2/2;215*2/2/2;30;50;20];
LB=[48*2/2/2;215*2/2/2;30;50;20];
%lb = [107;453;94;163;25];
%LB=[107;453;94;163;25];
%pause(5)
%[x,fval,exitflag] = ga(@TMD_nonlineardamploop,4,A,b,[],[],lb)
[x,fval,exitflag] = ga(@TMD_nonlineardamploop,5,A,b,[],[],LB,UB,[],[1,2,3,4,5])
function y = TMD_nonlineardamploop(x)
if ~isvector(x)
    error('Input must be a vector')
end
Kxt=x(1);
Kyt=x(2);
Kisayx=x(3)/1000;
Kisayy=x(4)/1000;
Kxt
Kyt
Kisayx
Kisayy
Ly=x(5)/100
m=30.;
Kx=1320; Ky=1320*90.0;
Kx=1320; Ky=1320*1.00;
wx=(Kx/m)^0.5;
%Kyt=13.5;
mt=1.5*2/2/2;
Lx=1.;
%Ly=1.;
Cx=2*0.025*4*(Kx/m)^0.5*m;
Cy=2*0.025*4*(Ky/m)^0.5*m;
%Kisayy=0.065;
%mt is for TMD
Cxt=2*Kisayx*(Kxt/mt)^0.5*mt; Cyt=2*Kisayy*(Kyt/mt)^0.5*mt;
%%stiffness & mass & damping
Beta=0.0;
%6063 7310 6160 6456
%309+15
%1270642
M=[m 0 0 0;0 m 0 0;0 0 mt 0;0 0 0 mt];
K=[Kx+Kxt 0 -Kxt 0;0 Ky+Kyt 0 -Kyt; -Kxt 0 Kxt 0; 0 -Kyt 0 Kyt];
C=[Cx+Cxt 0 -Cxt 0;0 Cy+Cyt 0 -Cyt; -Cxt 0 Cxt 0; 0 -Cyt 0 Cyt];
%eig(K*M^(-1))
wy=(Ky/m)^0.5;
%AAa=0.99;
www=1;
for F_wy=0.65*wx:0.1/2:1.25*wx
    for jjjj=1:2000.
        time=jjjj*0.02;
        
        AAa=0.65;
        
        %U_dotdotg(jjjj)=AAa*sin(1.*wy*time);
        U_dotdotg(jjjj)=AAa*sin(F_wy*time);
    end
    g=1;
    dt=0.02;
    
    
    %%---------------------------------------------
    %these parameters are for inital displacement, acceleration, velocity 
    x=zeros(length(U_dotdotg),4);
    v=zeros(length(U_dotdotg),4);
    a=zeros(length(U_dotdotg),4);
    %time 
    time=0:dt:dt*(length(U_dotdotg)-1);
    %these parameters are used for Newmark integration
    %aa = K*dt^2/2+C*dt;
    %bb = K*dt;
    %Keff=M+0.5*dt*C+dt^2/4*K;
    
    length(time);
    length(U_dotdotg);
    
    for i=1:length(U_dotdotg)-1
        
        i  ;
        %initial step calcultion
        disp(1,:)=x(i,:)+dt*v(i,:)+0.5*dt^2*a(i,:);
        velocity(1,:)=v(i,:)+dt*a(i,:);
        accel(1,:)=a(i,:);
        %initial guess of displacement, velocity, acceleration
        DD=disp(i,:)+dt*velocity(i,:)+dt^2/2*accel(i,:);
        VV=velocity(i,:)+dt*accel(i,:);
        AA=accel(i,:);
        criterion=10;
        j=1;
        %%ccc(1,:)=[0,0];
        Keff=M+0.5*dt*C+dt^2/4*K;
        eig(K^(-1)*M);
        while criterion>0.00005
            for ii=1:4
                xx=DD(3)-DD(1) ;
                yy=DD(4)-DD(2);
                VVx=VV(3)-VV(1);
                VVy=VV(4)-VV(2);
                fx1=Kxt*( Lx+xx  -(Lx^2+Lx*xx)/  ((Lx+xx   )^2+yy   ^2)^0.5)-Kxt*xx ; fx2=Kyt*(xx-   (Ly*xx)  /((Ly+yy)   ^2+xx^2   )^0.5) ;
                
                fxdamp1=Cxt*(VVx*(-yy^2/Lx^2+2*xx*yy^2/Lx^3)+VVy*(yy/Lx-xx*yy/Lx^2-yy^3/Lx^3+yy*xx^2/Lx^3)); 
                fxdamp2=Cyt*(VVx*(+xx^2/Ly^2-2*xx^2*yy/Ly^3)+VVy*(xx/Ly-xx*yy/Ly^2-xx^3/Lx^3+yy^2*xx/Ly^3));
                
                %%fXdamp1=Cxt*(VVx*(Lx+xx)^2+VVy*yy*(Lx+xx))/((Lx+xx)^2+yy^2)-Cxt*VVx;
                %%fXdamp2=Cyt*(VVy*(Ly+yy)*xx+VVx*xx^2)/((Ly+yy)^2+xx^2);
                %fx1=Kxt*( xx  +yy^2/2/Lx     -xx*yy^2/Lx^2)-Kxt*xx ; fx2=Kyt*(xx*yy/Ly     +xx^3/2/Ly^2      -xx*yy^2/Ly^2) ;
                
                fy1=Kxt*(yy   -(Lx*yy)  /((Lx+xx)^2+yy^2      )^0.5) ; fy2=  Kyt*( Ly+yy   -(Ly^2+Ly*yy) /((Ly+yy)   ^2+   xx^2)^0.5)-Kyt*yy ;
                
                fydamp1=Cxt*(VVy*(+yy^2/Lx^2-2*yy^2*xx/Lx^3)+VVx*(yy/Lx-yy*xx/Lx^2-yy^3/Ly^3+xx^2*yy/Lx^3));
                fydamp2=Cyt*(VVy*(-xx^2/Ly^2+2*yy*xx^2/Ly^3)+VVx*(xx/Ly-yy*xx/Ly^2-xx^3/Ly^3+xx*yy^2/Ly^3));
                
                %%fYdamp1=Cxt*(VVx*(Lx+xx)*yy+VVy*yy^2)/((Lx+xx)^2+yy^2);
                fYdamp2=Cyt*(VVy*(Ly+yy)^2+VVx*xx*(Ly+yy))/((Ly+yy)^2+xx^2)-Cyt*VVy;
                %fy1= Kxt*(xx*yy/Lx    +yy^3/2/Lx^2     -yy*xx^2/Lx^2) ;fy2=Kyt*( yy   +xx^2/2/Ly   -yy*xx^2/Ly^2)-Kyt*yy;
                if ii==1
                    F_N_S(ii)=-fx1-fx2+(-fxdamp1-fxdamp2);
                    Fe(i,ii)=-U_dotdotg(i)*g*cos(Beta);
                    Fe(i+1,ii)=-U_dotdotg(i+1)*g*cos(Beta);
                elseif ii==2
                    F_N_S(ii)=-fy1-fy2+(-fydamp1-fydamp2);
                    Fe(i,ii)=-U_dotdotg(i)*g*sin(Beta);
                    Fe(i+1,ii)=-U_dotdotg(i+1)*g*sin(Beta);
                elseif ii==3
                    F_N_S(ii)=fx1+fx2+(fxdamp1+fxdamp2);
                    Fe(i,ii)=-U_dotdotg(i)*g*cos(Beta);
                    Fe(i+1,ii)=-U_dotdotg(i+1)*g*cos(Beta);
                elseif ii==4
                    F_N_S(ii)=fy1+fy2+(fydamp1+fydamp2);
                    Fe(i,ii)=-U_dotdotg(i)*g*sin(Beta);
                    Fe(i+1,ii)=-U_dotdotg(i+1)*g*sin(Beta);    
                end 
            end
            %ccc(i+1,:)=[F_N_S(1)+Kxt*xx,F_N_S(2)+Kyt*yy];
            %ccc(i+1,:)=[Kxt*( Lx+xx  -(Lx^2+Lx*xx)/  ((Lx+xx   )^2+yy   ^2)^0.5),Kyt*(xx-   (Ly*xx)  /((Ly+yy)   ^2+xx^2   )^0.5)];
            %ccc(i+1,:)=[Kxt*(yy   -(Lx*yy)  /((Lx+xx)^2+yy^2      )^0.5),Kyt*( Ly+yy   -(Ly^2+Ly*yy) /((Ly+yy)   ^2+   xx^2)^0.5)];
            %%ccc(i+1,:)=[fydamp2,fYdamp2];
            Fe(i+1,:)=Fe(i+1,:)*M;
            
            deltaPhat=-AA*M +Fe(i+1,:) -VV*C -DD*K -F_N_S;
            
            deltaA=deltaPhat*Keff^(-1);
            DD=DD+dt^2*0.25*deltaA;
            VV=VV+dt*0.5*deltaA;
            AA=AA+deltaA;
            criterion=max(abs(deltaA));
            j=j+1;
            if j==95
                criterion=0.0;
            end
            NORM1(i)=max(abs(deltaPhat));
        end
        
        j    ;     
        iterationnum(i)=j;
        accel(i+1,:)=AA;
        velocity(i+1,:)=VV;
        disp(i+1,:)=DD;
        
    end
    
    MMM(www,:)=[max(abs(disp(1500:1999,1))),max(abs(disp(1500:1999,2)))];
    C_MMM(www)=www;
    C_MMM1(www)=F_wy;
    
    www=www+1;
    if mod(www,10)==0
        www;
    end
    %Axt=max(abs(disp(1500:1999,3)-disp(1500:1999,1)));
    %Ayt=max(abs(disp(1500:1999,4)-disp(1500:1999,2)));
    % AMP(www,:)=[max(abs(disp(1500:1999,1))),max(abs(disp(1500:1999,2))) ,Axt ,Ayt];
end
%plot(C_MMM1,[MMM(:,1),MMM(:,2)]),grid
%load 'disp_linear.txt'
non_linansw=[C_MMM1',MMM(:,1),MMM(:,2)];
%save('non_linansw_1.txt','non_linansw','-ASCII');
%load 'linansw.txt'
%plot(linansw(:,1)/6.63,[linansw(:,2),non_linansw(:,2),linansw(:,3),non_linansw(:,3)]),grid
%plot(AMP(:,1))
%plot(disp(2000:2500,3)-disp(2000:2500,1))
%plot(linansw(:,1)/6.63,[non_linansw(:,2),non_linansw(:,3)]),grid
%plot(time(:),[disp(:,2)-disp(:,4)]),grid
%plot(time(:),[ccc(:,1),ccc(:,2)]),grid
%plot(time(:),[ccc(:,1)]),grid
y =max( max([non_linansw(:,2),non_linansw(:,3)]));
y
end

3 个评论
显示 1更早的评论隐藏 1更早的评论

Rik 2023-6-27

编辑：Rik 2023-6-27

You forgot to ask a question; you just dumped uncommented code on us with a "help me". Have a read here and here. It will greatly improve your chances of getting an answer.

I edited your code, moving the part I suspect to be the calling script to the top. I tried running your code, but it took more than the 55 second limit.

Doha Ali 2023-6-27