Index exceeds the number of array elements. Index must not exceed 101.

Please help me fix the following error: "Index exceeds the number of array elements. Index must not exceed 101." thank you very much!

%% Phase plane
clear all
clc

%% conditions
A=[-0.8 0.1;0.15 -0.15];
A_h=[-0.4 0.2;0.1 0.04];
A_d=[0.2 0.5;-0.3 0.05];
A_u=[0.1;-0.1];


A_c=[-76.0089 387.8927;75.3589 -387.9427];

stime=1.0;
endtime=30;
h=0.01;     
t=-stime:h:endtime;
N_0=stime/h;
N_1=endtime/h;
N=N_0+N_1;
%bound=110;




for i=1:N+1 
          
    
          
    if i<=N_0+1;
           
          x1(:,i)=0.05*sin((i-N_0-1)*h); 
          x2(:,i)=0.05*cos((i-N_0-1)*h); 
         
         
    else
     d_1=floor(1.5+0.5*sin(1.5*i*h));
     d_2=floor(2+cos(2.5*i*h));
           
  if i<=4
      
     x1(:,i)=x1(:,i-1)+h*(A_c(1,1)*x1(:,i-1)+A_c(1,2)*x2(:,i-1));
               
     x2(:,i)=x2(:,i-1)+h*(A_c(2,1)*x1(:,i-1)+A_c(2,2)*x2(:,i-1));

          
  else
       x1(:,i)=x1(:,i-1)+h*(A_c(1,1)*x1(:,i-1)+A_c(1,2)*x2(:,i-1)+A_h(1,1)*x1(:,i-d_1-1)+A_h(1,2)*x2(:,i-d_1-1)+(1/2)*(A_d(1,1)*(x1(i-1)+x1(i-d_2-1)+2*(x1(i-d_2)+x1(i-d_2+1)))+A_d(1,2)*(x2(i-1)+x2(i-d_2-1)+2*(x2(i-d_2)+x2(i-d_2+1)))));
               
       x2(:,i)=x2(:,i-1)+h*(A_c(2,1)*x1(:,i-1)+A_c(2,2)*x2(:,i-1)+A_h(2,1)*x1(:,i-d_1-1)+A_h(2,2)*x2(:,i-d_1-1)+(1/2)*(A_d(2,1)*(x1(i-1)+x1(i-d_2-1)+2*(x1(i-2)+x1(i-3)+x1(i-4)))+A_d(2,2)*(x2(i-1)+x2(i-d_2-1)+2*(x2(i-2)+x2(i-3)+x2(i-4)))));
                
      

  end           
    end
  
end

%% Plotting Graphs

figure
grid on 
hold on
box on
plot(t,x1(1,:),'-g',t,x2(1,:),'r-','linewidth',2.5)
hold on
xlabel('Time (t)')
%ylabel('State responses')
legend('x_{1}(t)','x_{2}(t)');
xlim([0,30]);
ylim([-100,100]);

%

[Umar]

@Le,

To address your query regarding, “Please help me fix the following error: "Index exceeds the number of array elements. Index must not exceed 101." thank you very much!”, I have updated the code and was able to achieve your plot. Here is your updated code,

A = [-0.8 0.1; 0.15 -0.15];
A_h = [-0.4 0.2; 0.1 0.04];
A_d = [0.2 0.5; -0.3 0.05];
A_u = [0.1; -0.1];
A_c = [-76.0089 387.8927; 75.3589 -387.9427];
stime = 1.0;
endtime = 30;
h = 0.01;
t = -stime:h:endtime;
N_0 = stime / h;
N_1 = endtime / h;
N = N_0 + N_1;
x1 = zeros(2, N + 1); % Initialize x1 array
x2 = zeros(2, N + 1); % Initialize x2 array
for i = 1:N + 1
    if i <= N_0 + 1
        x1(:, i) = 0.05 * sin((i - N_0 - 1) * h);
        x2(:, i) = 0.05 * cos((i - N_0 - 1) * h);
    else
        d_1 = floor(1.5 + 0.5 * sin(1.5 * i * h));
        d_2 = floor(2 + cos(2.5 * i * h));
        
        if i <= 4
            x1(:, i) = x1(:, i - 1) + h * (A_c(1, 1) * x1(:, i - 1) + A_c(1, 2) * x2(:, i - 1));
            x2(:, i) = x2(:, i - 1) + h * (A_c(2, 1) * x1(:, i - 1) + A_c(2, 2) * x2(:, i - 1));
        else
            x1(:, i) = x1(:, i - 1) + h * (A_c(1, 1) * x1(:, i - 1) + A_c(1, 2) * x2(:, i - 1) + A_h(1, 1) * x1(:, i - d_1 - 1) + A_h(1, 2) * x2(:, i - d_1 - 1) + (1/2) * (A_d(1, 1) * (x1(i - 1) + x1(i - d_2 - 1) + 2 * (x1(i - d_2) + x1(i - d_2 + 1))) + A_d(1, 2) * (x2(i - 1) + x2(i - d_2 - 1) + 2 * (x2(i - d_2) + x2(i - d_2 + 1)))));
            x2(:, i) = x2(:, i - 1) + h * (A_c(2, 1) * x1(:, i - 1) + A_c(2, 2) * x2(:, i - 1) + A_h(2, 1) * x1(:, i - d_1 - 1) + A_h(2, 2) * x2(:, i - d_1 - 1) + (1/2) * (A_d(2, 1) * (x1(i - 1) + x1(i - d_2 - 1) + 2 * (x1(i - 2) + x1(i - 3) + x1(i - 4))) + A_d(2, 2) * (x2(i - 1) + x2(i - d_2 - 1) + 2 * (x2(i - 2) + x2(i - 3) + x2(i - 4)))));
        end
    end
end
%% Plotting Graphs
figure
grid on
hold on
box on
plot(t, x1(1, :), '-g', t, x2(1, :), 'r-', 'linewidth', 2.5)
xlabel('Time (t)')
legend('x_{1}(t)', 'x_{2}(t)');
xlim([0, 30]);
ylim([-100, 100]);

Please let me know if you have any further questions.

[Le]

thank you very much!

[Le]

The above code comes from the differential equation below. Can you help me check if I wrote the code correctly? Thank you for your dedicated help

[Torsten]

Don't try to code this on your own, you will fail with absolute certainty.

Instead, use MATLAB's "ddesd" to solve this system of delay differential equations.

The integrals can be included by solving the two additional differential equations

dI1/dt = x1(t)-x1(t-d2(t))*(1-d2'(t)), I1(0) = integral_{s=-3}^{s=0} x1(s) ds
dI2/dt = x2(t)-x2(t-d2(t))*(1-d2'(t)), I2(0) = integral_{s=-3}^{s=0} x2(s) ds

Then

I1(t) = integral_{s=t-d2(t)}^{s=t} x1(s) ds
I2(t) = integral_{s=t-d2(t)}^{s=t} x2(s) ds

in the course of the computation.

[Umar]

Hi @Le,

I do agree with @Torsten’s comments about, “ use MATLAB's "ddesd" to solve this system of delay differential equations”.

Brief summary: The script deli.m sets up a system of delay differential equations with defined functions for system equations, initial conditions, delays, parameters, and solver options. By executing the script and calling the ddesd solver, the delay differential equations are solved numerically. The solution is then extracted and plotted to visualize the behavior of the system over time.

function dxdt = ddex1de(t, x, Z, A, Ah, Ad)

    % Define the system of delay differential equations

    d1 = 1.5 + 0.5 * sin(1.5*t);

    d2 = 2 + cos(2*t);

    x1 = Z(:,1);

    x2 = Z(:,2);

    % trapz function

    dxdt = [A*x1 + Ah*x1 + Ad*trapz(t, x1) - x1; A*x2 + Ah*x2 + Ad*trapz(t, x2) - x2];

end

function Z = ddex1hist(t)

    % Define the initial conditions for x1 and x2

    Z = [0.05*sin(t); 0.05*cos(t)];

end

function d = ddex1delays(t)

    % Define the delay functions d1 and d2

    d = [1.5 + 0.5*sin(1.5*t); 2 + cos(2*t)];

end

% Define the parameters

A = 1.0;

Ah = 0.5;

Ad = 0.2;

% Define the options for the solver

options = ddeset('RelTol',1e-6,'AbsTol',1e-6);

% Solve the delay differential equations sol = ddesd(@ddex1de, @ddex1delays, @ddex1hist, [0, 5], options);

% Extract and plot the solution t = linspace(0, 5, 1000); y = deval(sol, t); plot(t, y(1,:), t, y(2,:)); legend('x1(t)', 'x2(t)'); xlabel('Time'); ylabel('Value'); title('Solution of Delay Differential Equations');

Please see attached plots and script.

[Torsten]

@Umar

dxdt = [A*x1 + Ah*x1 + Ad*trapz(t, x1) - x1; A*x2 + Ah*x2 + Ad*trapz(t, x2) - x2];

This cannot be used since t, x1 and x2 are scalar values so that trapz does not cover the interval [t-d2(t),t] for x1 and x2.

[Umar]

@Torsten,

Understood and completely aware that t, x1, and x2 should be defined as vectors representing the intervals over which you want to integrate. So, by providing arrays instead of scalars, you can utilize the trapz function effectively to compute the integral over the specified intervals for x1 and x2. But my point was trying to use MATLAB's "ddesd" function since never used it before and wanted to experiment with it. Also, I am aware that the trapz function expects vector inputs to compute the integral over a range, which is not feasible with scalar values.

[Le]

Thank @Umar and @Torsten. Do you think I wrote the code below correctly?

syms s t x1(s) x2(s)
x3(t) = int(x1,s,t-2-cos(2*t),t);
x4(t) = int(x2,s,t-2-cos(2*t),t);
diff(x3,t)

ans(t) = 

diff(x4,t)

ans(t) = 

int(0.05*sin(s),s,t-2-cos(2*t),t)

ans = 

int(0.05*cos(s),s,t-2-cos(2*t),t)

ans = 

sol = ddesd(@ddex1de,@ddex1delays,@ddex1hist,[0,30]);
plot(sol.x,[sol.y(1,:);sol.y(2,:)])

%xlim([0,30]);
%ylim([-0.5,0.5]);

function d = ddex1delays(t,y)
d = [t-1.5-0.5*sin(1.5*t);...
t-2-cos(2*t)];
end

function h = ddex1hist(t)
h = [0.05*sin(t);...
0.05*cos(t);...
(cos(cos(2*t)-t+2)-cos(t))/30;...
(sin(cos(2*t)-t+2)+sin(t))/30];
end

function dydt = ddex1de(t,y,Z)
dydt = [-0.8*y(1)+0.1*y(2)-0.4*Z(1,1)+0.2*Z(2,1)+0.2*y(3)+0.5*y(4);...
0.15*y(1)-0.15*y(2)+0.1*Z(1,1)+0.04*Z(2,1)-0.3*y(3)+0.05*y(4);...
y(1)-Z(1,2)*(2*sin(2*t)+1);...
y(2)-Z(2,2)*(2*sin(2*t)+1)];
end

[Torsten]

Seems to be the code I once tested here with MATLAB online.

The only difference is

h = [0.05*sin(t);...
0.05*cos(t);...
(cos(cos(2*t)-t+2)-cos(t))/30;...
(sin(cos(2*t)-t+2)+sin(t))/30];

I got

h = [0.05*sin(t);...
0.05*cos(t);...
(cos(cos(2*t)-t+2)-cos(t))/20;...
(sin(cos(2*t)-t+2)+sin(t))/20];

[Le]

Yes! thank you very much! @Torsten

[Umar]

Hi @Le,

Without alterations to code, I received the following results. In your code snippet, you have three functions related to Delay Differential Equations (DDEs) which are essential for solving DDE problems. From analysis of code, the ddex1delays function defines the delays in the DDE system. It returns a vector d representing the delays at time t with state variables y, The history function provides initial conditions or past values for the DDE system. It influences the system's behavior by incorporating past information into the current state and the last one ddex1de function defines the DDE system's differential equations. It computes the derivatives of the state variables y at time t with the history matrix Z. In nutshell, enabling the simulation and analysis of systems with delayed feedback. So, @Torsten is expert in field of dynamics and I would like his feedback if he can elaborate little bit about your recent comments, “The above code comes from the differential equation below. Can you help me check if I wrote the code correctly?”, and so far I haven’t see his feedback on it. However, my response to your comments “The above code comes from the differential equation below. Can you help me check if I wrote the code correctly?” is mentioned below.

Delay Functions:

The delay functions in the code are defined as:

        * d1(t) = t - 1.5 - 0.5sin(1.5t)
        * d2(t) = t - 2 - cos(2*t)

These functions correspond to the model's delay functions:

        * d11(t) = 1.5 + 0.5sin(1.5t)
        * d2(t) = 2 + cos(2*t)
        * My conclusion: The code implementation matches the specified delay functions.

Initial Conditions:

The initial conditions in the model are:

        * x1(t) = 0.05*sin(t)
        * x2(t) = 0.05*cos(t)

My conclusion: The code does not explicitly define these initial conditions but uses them in the history function ddex1hist(t). The history function ddex1hist(t) in the code corresponds to the initial conditions specified in the model.

Differential Equations:

The differential equations in the code are defined by the function ddex1de(t, y, Z). These equations involve the system dynamics and interactions based on the delayed states.

My conclusion: The equations in the code should be consistent with the mathematical model's system dynamics, considering the delays and interactions between variables. Hope, this helps.

[Torsten]

@Umar

Did you include the wrong graphics ? According to your code, you should get the same graphics as I do.

[Umar]

Hi @Torsten,

Honestly, I don’t think 🤔 so.

[Torsten]

So you didn't use this code to produce the graphics ?

sol = ddesd(@ddex1de,@ddex1delays,@ddex1hist,[0,30]);
plot(sol.x,[sol.y(1,:);sol.y(2,:)])

function d = ddex1delays(t,y)
d = [t-1.5-0.5*sin(1.5*t);...
t-2-cos(2*t)];
end
function h = ddex1hist(t)
h = [0.05*sin(t);...
0.05*cos(t);...
(cos(cos(2*t)-t+2)-cos(t))/30;...
(sin(cos(2*t)-t+2)+sin(t))/30];
end
function dydt = ddex1de(t,y,Z)
dydt = [-0.8*y(1)+0.1*y(2)-0.4*Z(1,1)+0.2*Z(2,1)+0.2*y(3)+0.5*y(4);...
0.15*y(1)-0.15*y(2)+0.1*Z(1,1)+0.04*Z(2,1)-0.3*y(3)+0.05*y(4);...
y(1)-Z(1,2)*(2*sin(2*t)+1);...
y(2)-Z(2,2)*(2*sin(2*t)+1)];
end

Obviously, code and graphics in your contribution don't fit together.

[Umar]

Hi @ Torsten,

Please feel more comfortable to glance through my screenshots as shown below. Hope, you will find something that I did not catch.

https://www.mathworks.com/matlabcentral/answers/uploaded_files/1746136/IMG_7305.jpeg

https://www.mathworks.com/matlabcentral/answers/uploaded_files/1745756/IMG_7291.jpeg

They are shared in this post.

Answer 1

[Star Strider]

I am not certain how to get these into a form that the MATLAB ODE integrators could use, so I did the next best thing, and calculated the derivatives and then integrated them numerically.

Try this —

%% conditions
A=[-0.8 0.1;0.15 -0.15];
A_h=[-0.4 0.2;0.1 0.04];
A_d=[0.2 0.5;-0.3 0.05];
A_u=[0.1;-0.1];

d_1 = @(t) floor(1.5+0.5*sin(1.5*t));
d_2 = @(t) floor(2+cos(2.5*t));



x1 = @(t) 0.05*sin(t); 
x2 = @(t) 0.05*cos(t);

x1_int = @(t,d) integral(@(t) x1(t), t-d_1(t), t, 'ArrayValued',1);
x2_int = @(t,d) integral(@(t) x2(t), t-d_2(t), t, 'ArrayValued',1);

tv = linspace(-3, 0, 250);

for k = 1:numel(tv)
    t = tv(k);
    xdot(:,k) = A*[x1(t); x2(t)] + A_h*[x1(t-d_1(t)); x2(t-d_2(t))] + A_d*[x1_int(t,d_1(t)); x2_int(t,d_2(t))];
end

xdot

xdot = 2x250

-0.0357 -0.0357 -0.0357 -0.0356 -0.0356 -0.0355 -0.0354 -0.0354 -0.0353 -0.0352 -0.0351 -0.0350 -0.0349 -0.0348 -0.0347 -0.0346 -0.0345 -0.0344 -0.0343 -0.0341 -0.0340 -0.0339 -0.0337 -0.0336 -0.0334 -0.0333 -0.0331 -0.0329 -0.0327 -0.0326 0.0029 0.0029 0.0028 0.0028 0.0028 0.0027 0.0027 0.0026 0.0026 0.0025 0.0025 0.0024 0.0024 0.0024 0.0023 0.0023 0.0022 0.0022 0.0021 0.0021 0.0020 0.0020 0.0019 0.0019 0.0018 0.0018 0.0017 0.0017 0.0016 0.0016

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

xm = cumtrapz(tv, xdot, 2)

xm = 2x250

0 -0.0004 -0.0009 -0.0013 -0.0017 -0.0021 -0.0026 -0.0030 -0.0034 -0.0039 -0.0043 -0.0047 -0.0051 -0.0055 -0.0060 -0.0064 -0.0068 -0.0072 -0.0076 -0.0080 -0.0084 -0.0089 -0.0093 -0.0097 -0.0101 -0.0105 -0.0109 -0.0113 -0.0117 -0.0121 0 0.0000 0.0001 0.0001 0.0001 0.0002 0.0002 0.0002 0.0003 0.0003 0.0003 0.0004 0.0004 0.0004 0.0004 0.0005 0.0005 0.0005 0.0005 0.0006 0.0006 0.0006 0.0006 0.0007 0.0007 0.0007 0.0007 0.0008 0.0008 0.0008

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

figure
yyaxis left
plot(tv, xdot(1,:))
ylabel('$\dot{x_1}(t)$', 'Interpreter','LaTeX')
yyaxis right
plot(tv, xdot(2,:))
ylabel('$\dot{x_2}(t)$', 'Interpreter','LaTeX')
grid
xlabel('Time')

figure
yyaxis left
plot(tv, xm(1,:))
ylabel('x_1(t)')
yyaxis right
plot(tv, xm(2,:))
ylabel('x_2(t)')
grid
xlabel('Time')

The derivatives look a bit strange to me, however you have access to the paper and know what it wants, so I defer to you for any necessary corrections.

I do not see anywhere in the excerpt from the paper where ‘A_c’ appears, so I did not use it.

.

[Le]

Thank yor very much! How do the orbitals x_1(t) and x_2(t) appear on the interval [0,30]?

[Star Strider]

The paper excerpt defined t as going from to 0 so I used those in my analysis.

Changing them to is striaghtforward —

%% conditions
A=[-0.8 0.1;0.15 -0.15];
A_h=[-0.4 0.2;0.1 0.04];
A_d=[0.2 0.5;-0.3 0.05];
A_u=[0.1;-0.1];

d_1 = @(t) floor(1.5+0.5*sin(1.5*t));
d_2 = @(t) floor(2+cos(2.5*t));



x1 = @(t) 0.05*sin(t); 
x2 = @(t) 0.05*cos(t);

x1_int = @(t,d) integral(@(t) x1(t), t-d_1(t), t, 'ArrayValued',1);
x2_int = @(t,d) integral(@(t) x2(t), t-d_2(t), t, 'ArrayValued',1);

tv = linspace(0, 30, 3000);

for k = 1:numel(tv)
    t = tv(k);
    xdot(:,k) = A*[x1(t); x2(t)] + A_h*[x1(t-d_1(t)); x2(t-d_2(t))] + A_d*[x1_int(t,d_1(t)); x2_int(t,d_2(t))];
end

xdot

xdot = 2x3000

0.0109 0.0358 0.0358 0.0359 0.0359 0.0359 0.0359 0.0359 0.0359 0.0358 0.0358 0.0358 0.0358 0.0358 0.0357 0.0357 0.0357 0.0356 0.0356 0.0356 0.0355 0.0355 0.0354 0.0353 0.0353 0.0352 0.0351 0.0351 0.0350 0.0349 -0.0064 -0.0033 -0.0033 -0.0033 -0.0033 -0.0032 -0.0032 -0.0032 -0.0031 -0.0031 -0.0031 -0.0030 -0.0030 -0.0030 -0.0029 -0.0029 -0.0029 -0.0028 -0.0028 -0.0028 -0.0027 -0.0027 -0.0026 -0.0026 -0.0026 -0.0025 -0.0025 -0.0025 -0.0024 -0.0024

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

xm = cumtrapz(tv, xdot, 2)

xm = 2x3000

0 0.0002 0.0006 0.0010 0.0013 0.0017 0.0020 0.0024 0.0027 0.0031 0.0035 0.0038 0.0042 0.0045 0.0049 0.0053 0.0056 0.0060 0.0063 0.0067 0.0070 0.0074 0.0077 0.0081 0.0084 0.0088 0.0092 0.0095 0.0099 0.0102 0 -0.0000 -0.0001 -0.0001 -0.0001 -0.0002 -0.0002 -0.0002 -0.0003 -0.0003 -0.0003 -0.0004 -0.0004 -0.0004 -0.0005 -0.0005 -0.0005 -0.0005 -0.0006 -0.0006 -0.0006 -0.0007 -0.0007 -0.0007 -0.0007 -0.0008 -0.0008 -0.0008 -0.0008 -0.0009

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

figure
yyaxis left
plot(tv, xdot(1,:))
ylabel('$\dot{x_1}(t)$', 'Interpreter','LaTeX')
yyaxis right
plot(tv, xdot(2,:))
ylabel('$\dot{x_2}(t)$', 'Interpreter','LaTeX')
grid
xlabel('Time')
title('Derivatives')

figure
yyaxis left
plot(tv, xm(1,:))
ylabel('x_1(t)')
yyaxis right
plot(tv, xm(2,:))
ylabel('x_2(t)')
grid
xlabel('Time')
title('Integrated Result')

Note that I used yyaxis so these plotted results will appear different from the others.

.

[Le]

Ok, thank you very much!

[Star Strider]

My pleasure!

[Le]

Do you have any code for the event-triggered control @Star Strider?

[Star Strider]

This version of your code does not perrmit events, although it possibly could if I knew what you want to do.

[Le]

I want to write control code for the system class (1) below but can't? do you help me? I am grateful to you

[Umar]

@Le,

This is a very complex and interesting model in field of dynamics and not easy to implement. I am aware if you integrate the state feedback controller and the event-triggered scheme into the Roesser model, a comprehensive control strategy can be developed to regulate the behavior of the 2-D system with delays effectively which enhances the system's stability and performance, making it a crucial aspect of control system design and implementation.

[Umar]

Hi @Le,

To address your query regarding, “I want to write control code for the system class (1) below but can't? do you help me? I am grateful to you”

I am providing code framework for the control code for the system class based on the Roesser model. It defines the system parameters, state feedback controller, event-triggered scheme parameters, initial conditions, time vector, and implements the control system using the Roesser model with delays. The results are then plotted to visualize the behavior of the system over time. Here is the framework of code,

% Define system parameters with random perturbations

A = [0.01 0.504; -0.5 0.14] + 0.1*randn(2);

B = [0; 0.25] + 0.05*randn(2);

C = [1 0];

D = 0;

delay = 1; % Delay value

% State feedback controller with random perturbations

K = [2.393 -1.472] + 0.1*randn(1,2);

% Event-triggered scheme parameters with random perturbations

alpha_u = 0.5 + 0.1*randn;

alpha_v = 0.8 + 0.1*randn;

Ph = 3 + 0.1*randn;

Pv = 5 + 0.1*randn;

% Initial conditions with random perturbations

x1_0 = 0.5 + 0.1*randn; % Updated random value for x1_0

x2_0 = 0 + 0.1*randn;

xv_0 = 0.3 + 0.1*randn; % Updated random value for xv_0

% Time vector

t = 0:0.01:10; % Define the time span

% Initialize arrays to store state values

x1 = zeros(size(t));

x2 = zeros(size(t));

xv = zeros(size(t));

u = zeros(size(t));

% Implement the control system with random perturbations

for i = 2:length(t)

    % Calculate error terms

    e1 = x1(i-1) - x1_0;

    e2 = x2(i-1) - xv_0;

    % Calculate the control input based on the event-triggered scheme

    u(i) = K * [x1(i-1); xv(i-1)] + 0.1*randn; % Updated random value for control input u

    % Update states using the Roesser model with delays

    x1(i) = A(1,1)*x1(i-1) + A(1,2)*x2(i-1) + B(1)*u(i-delay);

    x2(i) = A(2,1)*x1(i-1) + A(2,2)*x2(i-1) + B(2)*u(i-delay);

    xv(i) = x1(i-delay);

end

% Plot the results with random perturbations

figure;

subplot(2,1,1);

plot(t, x1, 'b', t, xv, 'r');

xlabel('Time');

ylabel('x1, xv');

legend('x1', 'xv');

title('State Variables x1 and xv over Time with Random Perturbations');

subplot(2,1,2);

plot(t, u, 'g');

xlabel('Time');

ylabel('Control Input u');

title('Control Input u over Time with Random Perturbations');

Please see attached plots.

Hope this answers your question. Please customize the code based on your preferences.

[Le]

Thank you for your enthusiasm in helping me @Umar. How to get the 2D orbits of x^h(i,j), x^v(i,j) and the horizontal event-triggered instants and corresponding event-triggered interval? Figures as I attached the file below.

[Umar]

Hi @Le,

Generate 2D orbits for x^h(i,j) and x^v(i,j) by defining the range of i and j values using linspace to create a set of points, create a grid of points (i,j) using meshgrid, calculate the values of x^h(i,j) and x^v(i,j) based on the given functions.Finally, plot the 2D orbits of x^h(i,j) and x^v(i,j) using surf to visualize the orbits.To determine the horizontal event-triggered instants and their intervals.Here is an example to help you out, First initialize parameters such as the number of iterations, step size, and threshold for event triggering, then create Arrays to store the horizontal and vertical orbits, event-triggered instants, and intervals. Then, use iterative calculations to generate 2D orbits, event triggering is checked based on the threshold, and relevant data is stored.

[Le]

This problem is relatively difficult for me. Can you write me a sample code? Thank you with all my heart!

[Umar]

Hi @Le,

How about you implement the framework and then I will take a look into it to help you further.

[Le]

Hi @Umar : Below is a code about state feedback controller.

clc
clear all;
close all;

%Closed-loop system 
A=[0.5 1;0.45 0.85];
A_d=[0.04 0.1;0.1 0.2];
B=[1.5;1.3];
D=[0.5;0.4];
g=7;
K=[-0.3202 -0.7155];


A_c=A+B*K;

n=40;
v=zeros(1,n);
w=zeros(1,n);
Xh=zeros(n,n); Xv=zeros(n,n);U=zeros(n,n);W=zeros(n,n);

for i=1:1:10
    Xh(1,i)=0;
end
for i=11:1:n-1
    Xh(1,i)=0;
  
end
for i=1:1:10
    Xv(i,1)=0;
end
for i=11:1:n-1
    Xv(i,1)=0; 
end
x1=zeros(2,1);x0=zeros(2,1);xt=zeros(2,1);xd=ones(2,1);
for i=1:1:n
    v(1,i)=i;
    w(1,i)=i;
end

for s=1:1:4
    for i=s:1:4
        x0(1,1)=Xh(s,i);x0(2,1)=Xv(s,i);
        x1=(A_c*x0)+(D*(9*exp(-0.065*pi*(s+i))));
        Xh(s+1,i)=x1(1,1); 
        Xv(s,i+1)=x1(2,1);
        U(s,i)=K*x0;
        W(s,i)=9*exp(-0.065*pi*(s+i));
    end
end
for s=1:1:4
    for i=5:1:n-1
    x0(1,1)=Xh(s,i);x0(2,1)=Xv(s,i);
    xt(1,1)=Xh(s,i-3); xt(2,1)=Xv(s,i-3);
    x1=(A_c*x0)+(A_d*xt)+(D*(9*exp(-0.065*pi*(s+i))));    
    Xh(s+1,i)=x1(1,1); 
    Xv(s,i+1)=x1(2,1);
    U(s,i)=K*x0;
    end
end

for s=5:1:n-1
    for i=s:1:n-1
     x0(1,1)=Xh(s,i);
     x0(2,1)=Xv(s,i);
    xt(1,1)=Xh(s-3,i-3);
    xt(2,1)=Xv(s-3,i-3);
    x1=(A_c*x0)+(A_d*xt)+(D*(9*exp(-0.065*pi*(s+i))));    
    Xh(s+1,i)=x1(1,1); 
    Xv(s,i+1)=x1(2,1);
    U(s,i)=K*x0;
    end
end

for i=1:1:4
    for s=i:1:4
        x0(1,1)=Xh(s,i);x0(2,1)=Xv(s,i);
        x1=(A_c*x0)+(D*(9*exp(-0.065*pi*(s+i))));
        Xh(s+1,i)=x1(1,1); 
        Xv(s,i+1)=x1(2,1);
        U(s,i)=K*x0;
    end
end

for i=1:1:4
    for s=5:1:n-1
    x0(1,1)=Xh(s,i);
    x0(2,1)=Xv(s,i);
    xt(1,1)=Xh(s-3,i);
    xt(2,1)=Xv(s-3,i);
    x1=(A_c*x0)+(A_d*xt)+(D*(9*exp(-0.065*pi*(s+i))));    
    Xh(s+1,i)=x1(1,1); 
    Xv(s,i+1)=x1(2,1);
    U(s,i)=K*x0;
    end
end
for i=5:1:n-1
    for s=i:1:n-1
     x0(1,1)=Xh(s,i);
     x0(2,1)=Xv(s,i);
    xt(1,1)=Xh(s-3,i-3);
    xt(2,1)=Xv(s-3,i-3);
    x1=(A_c*x0)+(A_d*xt)+(D*(9*exp(-0.065*pi*(s+i))));    
    Xh(s+1,i)=x1(1,1);
    Xv(s,i+1)=x1(2,1);
    U(s,i)=K*x0;
    end
end




figure(1)
surf(v,w,Xh)
xlabel('i');
ylabel('j');
zlabel('x^h(i,j)');
axis([0 40 0 40 -1 3])
mycolors = [0 1 1; 1 0 0; 1 0 0];
colormap(mycolors);
%set(gca,'FontSize',18);

figure(2)
surf(v,w,Xv)
xlabel('i');
ylabel('j');
zlabel('x^v(i,j)');
axis([0 40 0 40 -1 3])
mycolors = [1 1 0; 1 0 0; 1 0 0];
colormap(mycolors);


figure(3)
surf(v,w,U)
xlabel('i');
ylabel('j');
zlabel('u(i,j)');
axis([0 40 0 40 -3 1])
mycolors = [0 1 1; 1 1 0; 1 0 1];
colormap(mycolors);

Can you help me replace state feedback controller by event trigggered control u(i,j) =K[x^h(i_p,j);x^v(i,j_q)]?

Thank you very much!

[Umar]

Hi @Le,

To adress your query regarding, “can you help me replace state feedback controller by event trigggered control u(i,j) =K[x^h(i_p,j);x^v(i,j_q)]?”

To replace the state feedback controller with event-triggered control in your given code, you need to make several modifications: first define the triggering conditions by assuming that the triggering condition is based on the difference between the current state and the previous state. Then, modify the control update only when the triggering condition is satisfied. This can be done by adding an if statement to check if the triggering condition is met. If it is, we will update the control action; otherwise, keep the previous control action. Since event-triggered control is being utilized, the state variables will be updated only when the control action is updated. Therefore, modifying the code to update the state variables accordingly. Here is your modified code:

% Closed-loop system

A = [0.5 1; 0.45 0.85];

A_d = [0.04 0.1; 0.1 0.2];

B = [1.5; 1.3];

D = [0.5; 0.4];

g = 7;

K = [-0.3202 -0.7155];

A_c = A + B * K;

n = 40;

v = zeros(1, n);

w = zeros(1, n);

Xh = zeros(n, n);

Xv = zeros(n, n);

U = zeros(n, n);

W = zeros(n, n);

prev_control = zeros(2, 1); % Initialize previous control action

for i = 1:1:10

    Xh(1, i) = 0;

end

for i = 11:1:n-1

    Xh(1, i) = 0;

end

for i = 1:1:10

    Xv(i, 1) = 0;

end

for i = 11:1:n-1

    Xv(i, 1) = 0;

end

x1 = zeros(2, 1);

x0 = zeros(2, 1);

xt = zeros(2, 1);

xd = ones(2, 1);

for i = 1:1:n

    v(1, i) = i;

    w(1, i) = i;

end

for s = 1:1:4

    for i = s:1:4

        x0(1, 1) = Xh(s, i);

        x0(2, 1) = Xv(s, i);

        x1 = (A_c * x0) + (D * (9 * exp(-0.065 * pi * (s + i))));

        Xh(s + 1, i) = x1(1, 1);

        Xv(s, i + 1) = x1(2, 1);

        % Check triggering condition and update control action

        if norm(x1 - prev_control) > threshold

            U(s, i) = K * x0;

            prev_control = U(s, i);

        else

            U(s, i) = prev_control;

        end

        W(s, i) = 9 * exp(-0.065 * pi * (s + i));

    end

end

After the last end statement, the rest of your code remains the same. So, in this modified code, you can see an if statement is added to check the triggering condition before updating the control action. The triggering condition is based on the difference between the current state x1 and the previous control action prev_control. If the difference exceeds a certain threshold, control action is updated, otherwise, you need to keep the previous control action. Please note that you need to define the threshold value based on your specific requirements. The threshold determines how much the state needs to change before triggering a control update. By implementing this event-triggered control, you can reduce the computational burden and improve the efficiency of the control system. The control action will only be updated when necessary, leading to better resource utilization and potentially improved system performance. I hope this helps! Let me know if you have any further questions.

Answer 2

[Voss]

The terms x1(i-d_2+1) and x2(i-d_2+1) will cause that error because d_2 is 1.

[Le]

CAN YOU RE-EDIT THE ABOVE CODE SO THAT IT CAN RUN! THANK YOU SO MUCH

[Voss]

Sure, I can edit it to run in a million different ways. Which way would be the right way? I have no idea.

[Le]

You can use the method of keeping the d_2=floor(2+ (cos(2.5*i*h))) fixed and adjusting the

x1(:,i)=x1(:,i-1)+h*(A_c(1,1)*x1(:,i-1)+A_c(1,2)*x2(:,i-1)+A_h(1,1)*x1(:,i-d_1-1)+A_h(1,2)*x2(:,i-d_1-1)+(1/2)*(A_d(1,1)*(x1(i-1)+x1(i-d_2-1)+A_d(1,2)*(x2(i-1)+x2(i-d_2-1)))));

and

x2(:,i)=x2(:,i-1)+h*(A_c(2,1)*x1(:,i-1)+A_c(2,2)*x2(:,i-1)+A_h(2,1)*x1(:,i-d_1-1)+A_h(2,2)*x2(:,i-d_1-1)+(1/2)*(A_d(2,1)*(x1(i-1)+x1(i-d_2-1))+A_d(2,2)*(x2(i-1)+x2(i-d_2-1))));

??

[Torsten]

You can use the method of keeping the d_2=floor(2+ (cos(2.5*i*h))) fixed and adjusting the ...

@Voss asked for the correct adjustment. It's only you who can tell us.