what is wrong in my code , there is a mistake in it because it is not converged

Question

Doha Ali 2023-5-1

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1955754-what-is-wrong-in-my-code-there-is-a-mistake-in-it-because-it-is-not-converged

编辑： Doha Ali 2023-5-5

CHICHI0968.mat

clear;clc;warning off;
global Kx Ky Kxt Kyt m mt Cx Cy Cxt Cyt xg Lx Ly j Beta  wx w
Kx=1500;Ky=1500;Kxt=150;Kyt=150;m=10;mt=1;
Lx=1.;
Ly=1.;
wx=(Kx/m)^0.5;
%mt is for TMD
Cx=2*0.02*(Kx/m)^0.5*m; Cy=2*0.02*(Ky/m)^0.5*m; Cxt=2*0.02*(Kxt/mt)^0.5*mt; Cyt=2*0.02*(Kyt/mt)^0.5*mt; 
Beta=pi/6;
w=1;
%--------------------------------------------------------------------------
A=.03;
dt=.01;tf=30.;t=0:dt:tf;n=tf/dt;%tsp=time step; tf=final time;
load('CHICHI0968.mat');%CHICHI0968.mat-max of elcentro=0.3487
%ug=9.81*CHICHI0968(2,1:n+1);%%ELCENTRO_NS0348;A*(wx*w)^2*sin(wx*w*t);
ug=A*(wx*w)^2*sin(wx*w*t);
xg=ug;
%--------------------------------------------------------------------------
x0N=[0 0 0 0 0 0 0 0];x0L=x0N;xjN(1,:)=x0N;xjL(1,:)=x0L;
for j=1:n
    %for j=1:3500   
    tint=dt*[j-1 j];%tint=time interval
    [tN,xN] = ode45(@nonlinearmodel,tint,x0N);
    xjN(j+1,:)=xN(length(tN),:);
    x0N=xN(length(tN),:);
    [tL,xL] = ode45(@linearmodel,tint,x0L);
    xjL(j+1,:)=xL(length(tL),:);
    x0L=xL(length(tL),:);
    j;
end
%--------------------------------------------------------------------------
ux1=[xjN(:,1) xjL(:,1)];ux2=[xjN(:,2) xjL(:,2)];
uy1=[xjN(:,3) xjL(:,3)];uy2=[xjN(:,4) xjL(:,4)];
%ux1=[xjL(:,1)];ux2=[xjL(:,2)];
%uy1=[xjN(:,3) xjL(:,3)];uy2=[xjN(:,4) xjL(:,4)];
global Kx Ky Kxt Kyt m mt Cx Cy Cxt Cyt xg Lx Ly j Beta wx
% Define parameters
Kx=1500; Ky=1500; Kxt=150; Kyt=150; m=10; mt=1;
Lx=1.; Ly=1.;
wx=(Kx/m)^0.5;
Cx=2*0.02*(Kx/m)^0.5*m; Cy=2*0.02*(Ky/m)^0.5*m; Cxt=2*0.02*(Kxt/mt)^0.5*mt; Cyt=2*0.02*(Kyt/mt)^0.5*mt; 
Beta=pi/6;
w=1;
% Define optimization problem
fun = @(A) cost_function(A);
x0 = 0.03;
lb = 0.01;
ub = 0.1;
options = optimoptions('fmincon','Display','iter','Algorithm','sqp');
[A_opt, J_opt] = fmincon(fun, x0, [], [], [], [], lb, ub, [], options);
figure;
subplot(2,1,1);
plot(t, ux1(:,1), 'b', t, ux1(:,2), 'r');
xlabel('Time (s)');
ylabel('Displacement (m)');
legend('Nonlinear Model', 'Linear Model');
title('Displacement of Main Mass (ux1)');
subplot(2,1,2);
plot(t, uy1(:,1), 'b', t, uy1(:,2), 'r');
xlabel('Time (s)');
ylabel('Displacement (m)');
legend('Nonlinear Model', 'Linear Model');
title('Displacement of TMD (uy1)');
% Define function for cost function
function J = cost_function(A)
global Kx Ky Kxt Kyt m mt Cx Cy Cxt Cyt xg Lx Ly j Beta  wx w
dt=0.01;
tf=30.;
t=0:dt:tf;
n=tf/dt;
load('CHICHI0968.mat');
ug=A*(wx*w)^2*sin(wx*w*t);
xg=ug;
x0N=[0 0 0 0 0 0 0 0];
x0L=x0N;
xjN(1,:)=x0N;
xjL(1,:)=x0L;
for j=1:n   
    tint=dt*[j-1 j];
    [tN,xN] = ode45(@nonlinearmodel,tint,x0N);
    xjN(j+1,:)=xN(length(tN),:);
    x0N=xN(length(tN),:);
    [tL,xL] = ode45(@linearmodel,tint,x0L);
    xjL(j+1,:)=xL(length(tL),:);
    x0L=xL(length(tL),:);
end
ux1=[xjN(:,1) xjL(:,1)];
ux2=[xjN(:,2) xjL(:,2)];
uy1=[xjN(:,3) xjL(:,3)];
uy2=[xjN(:,4) xjL(:,4)];
% Define cost function as the maximum displacement of the main structure
J = max(abs(ux1(:,1)));
end
% Define differential equation for nonlinear model
function dx=nonlinearmodel(t,x)
global Kx Ky Kxt Kyt m mt Cx Cy Cxt Cyt xg Lx Ly j Beta 
dx = zeros(8,1);
dx(1)=x(5);
dx(2)=x(6);
dx(3)=x(7);
dx(4)=x(8);
dx(5)=1/m*(-Cx*x(5)-Kx*x(1)-Cxt*(x(5)-x(6))+Kxt*((x(2)-x(1))+(x(4)-x(3))^2/(2*Lx)-(x(2)-x(1))*(x(4)-x(3))^2/Lx^2))...
    +1/m *Kyt*((x(2)-x(1))*(x(4)-x(3))/Ly -(x(2)-x(1))*(x(4)-x(3))^2/Ly^2+(x(2)-x(1))^3/(2*Ly^2))-xg(j)*cos(Beta);
dx(6)=1/mt*(Cxt*(x(5)-x(6))-Kxt*((x(2)-x(1))+(x(4)-x(3))^2/(2*Lx)-(x(2)-x(1))*(x(4)-x(3))^2/Lx^2))-xg(j)*cos(Beta);
dx(7)=1/m*(-Cy*x(7)-Ky*x(3)-Cyt*(x(7)-x(8))+Kyt*((x(4)-x(3))+(x(2)-x(1))^2/(2*Ly)-(x(4)-x(3))*(x(2)-x(1))^2/Ly^2))...
    +1/m*Kxt*((x(4)-x(3))*(x(2)-x(1))/Lx -(x(4)-x(3))*(x(2)-x(1))^2/Lx^2+(x(4)-x(3))^3/(2*Lx^2))-xg(j)*sin(Beta);
dx(8)=1/mt*(Cxt*(x(7)-x(8))-Kyt*((x(4)-x(3))+(x(2)-x(1))^2/(2*Ly)-(x(4)-x(3))*(x(2)-x(1))^2/Ly^2))...
    -1/mt*Kxt*((x(4)-x(3))*(x(2)-x(1))/Lx -(x(4)-x(3))*(x(2)-x(1))^2/Lx^2+(x(4)-x(3))^3/(2*Lx^2))-xg(j)*sin(Beta);
end
% Define differential equation for linear model
function dx=linearmodel(t,x)
global Kx Ky Kxt Kyt m mt Cx Cy Cxt Cyt xg Lx Ly j Beta 
dx = zeros(8,1);
dx(1)=x(5);
dx(2)=x(6);
dx(3)=x(7);
dx(4)=x(8);
dx(5)=1/m*(-Cx*x(5)-Kx*x(1)-Cxt*(x(5)-x(6))+Kxt*((x(2)-x(1))))-xg(j)*cos(Beta);
dx(6)=1/mt*(Cxt*(x(5)-x(6))-Kxt*((x(2)-x(1))))-xg(j)*cos(Beta);
dx(7)=1/m*(-Cy*x(7)-Ky*x(3)-Cyt*(x(7)-x(8))+Kyt*((x(4)-x(3))))-xg(j)*sin(Beta);
dx(8)=1/mt*(Cxt*(x(7)-x(8))-Kxt*(x(4)-x(3)))-xg(j)*sin(Beta);
end