필터 지우기
필터 지우기

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

조회 수: 1 (최근 30일)
Doha Ali
Doha Ali 2023년 5월 1일
편집: Doha Ali 2023년 5월 5일
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
Torsten
Torsten 2023년 5월 5일
편집: Torsten 2023년 5월 5일
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.
Doha Ali
Doha Ali 2023년 5월 5일
편집: Doha Ali 2023년 5월 5일
this is the mat . files and the load , please help me

댓글을 달려면 로그인하십시오.

답변 (0개)

카테고리

Help CenterFile Exchange에서 Nonlinear Control에 대해 자세히 알아보기

제품


릴리스

R2009b

Community Treasure Hunt

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

Start Hunting!

Translated by