Clear Filters
Clear Filters

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

1 view (last 30 days)
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
  2 Comments
Torsten
Torsten on 5 May 2023
Edited: Torsten on 5 May 2023
It will be hard to find the reason for non-convergence even with the .mat-files you have to load, but without them, it will be impossible.

Sign in to comment.

Answers (0)

Categories

Find more on Nonlinear Control in Help Center and File Exchange

Products


Release

R2009b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!