not enough input arguments

hi, I'm trying to write a code about cstr simulation and my loop cant continue due to this error, I'm sure it has what it needs for just first iteration but it does not calculate after the k1 part, can you help me please.

clear; clc;
Ts = 0.1;
%t = [1:0.1:1000];
Da1 = 3;
Da2 = 0.5;
Da3 = 1;
x1(1) = 0.5;
x2(1) = 0.5;
x3(1) = 1.5;
u = 0.45;
d2 = 1;
for n=1:1000
    k1x1(n) = Ts*f1_func(x1(n),x2(n),Da1,Da2)
    k1x2(n) = Ts*f2_func(x1(n),x2(n),Da1,Da2,Da3,d2,u)
    k1x3(n) = Ts*f3_func(x2(n),x3(n),Da3,d2)
    k2x1(n) = Ts*f1_func(x1(n)+k1x1(n)/2,x2(n)+k1x2(n)/2,u);
    k2x2(n) = Ts*f2_func(x1(n)+k1x1(n)/2,x2(n)+k1x2(n)/2,u);
    k2x3(n) = Ts*f3_func(x1(n)+k1x1(n)/2,x2(n)+k1x2(n)/2,x3(n)+k1x3(n)/2);
    k3x1(n) = Ts*f1_func(x1(n)+k2x1(n)/2,x2(n)+k2x2(n)/2,u);
    k3x2(n) = Ts*f2_func(x1(n)+k2x1(n)/2,x2(n)+k2x2(n)/2,u);
    %k3x3
    k4x1(n) = Ts*f1_func(x1(n)+k3x1(n),x2(n)+k3x2(n),u);
    k4x2(n) = Ts*f2_func(x1(n)+k3x1(n),x2(n)+k3x2(n),u);
    %k4x3
    x1(n+1) = x1(n)+1/6*(k1x1(n)+2*k2x1(n)+2*k3x1(n)+k4x1(n));
    x2(n+1) = x2(n)+1/6*(k1x2(n)+2*k2x2(n)+2*k3x2(n)+k4x2(n));
    %x3
end

k1x1 = -0.0875

k1x2 = 0.1075

k1x3 = -0.1250

Not enough input arguments.

Error in f1_func (line 3)
dx1 = 1 - x1(t) - Da1*x1(t) + Da2*x2(t)*x2(t);

Answer 1

[Walter Roberson]

You declare f1_func as expecting f1_func(x1,x2,Da1,Da2) but starting from

    k2x1(n) = Ts*f1_func(x1(n)+k1x1(n)/2,x2(n)+k1x2(n)/2,u);

you only pass in three parameters.

Also

        k1x2(n) = Ts*f2_func(x1(n),x2(n),Da1,Da2,Da3,d2,u)

you are passing in 7 parameters there, but

        k2x2(n) = Ts*f2_func(x1(n)+k1x1(n)/2,x2(n)+k1x2(n)/2,u);

you are passing in 3 parameters there.

[Halil Ibrahim]

i think i understand now but one more thing. so when i define the first function with (x1,x2,Da1,Da2,u) these parameters and later in my script when i input (x1(n)+k1x1(n)/2,x2(n)+k1x2(n)/2,u), is it directly replacing (x1) with (x1(n)+k1x1(n)/2), x2 with (x2(n)+k1x2(n)/2) and so on ?

[Sam Chak]

Hi @Halil Ibrahim

If the dynamics (ordinary differential equations) of the Continuous Stirred-Tank Reactor (CSTR) are given, would you like to run the simulation using the built-in ode45() solver? This way, you can focus on ensuring that the CSTR equations and parameters are designed or described correctly.

The code snippet is as follows:

%% Simulation settings and call ode45 solver
tspan   = 0:0.1:1000;       % simulation time from 0 to 1000
x10     = 0.5;
x20     = 0.5;
x30     = 1.5;
x0      = [x10; x20; x30];  % initial values of states x1, x2, x3
[t, x]  = ode45(@cstr, tspan, x0);

%% Plot results
plot(t, x), grid, xlabel Time

%% Continuous Stirred-Tank Reactor
function dxdt = cstr(t, x)
    dxdt    = zeros(3, 1);
    
    % Definitions of the states
    cA      = x(1); % concentration of A
    cB      = x(2); % concentration of B
    cC      = x(3); % concentration of C
    
    % Constants
    fAi     = ...;  % molar flow rate inlet of species A
    fAo     = ...;  % molar flow rate outlet of species A
    vA      = ...;  % stoichiometric coefficient
    rA      = ...;  % reaction rate
    
    % Parameters that depend on the constants or states (if any)
    p1      = ...;
    p2      = ...;
    p3      = ...;
    
    % Differential equations
    dxdt(1) = ...;
    dxdt(2) = ...;
    dxdt(3) = ...;
end

[Walter Roberson]

when i input (x1(n)+k1x1(n)/2,x2(n)+k1x2(n)/2,u), is it directly replacing (x1) with (x1(n)+k1x1(n)/2), x2 with (x2(n)+k1x2(n)/2) and so on ?

Yes. The replacements are positional. Whatever is passed first will be associated with the first variable in the function() definition; whatever is passed second will be associated with the second variable in the function() definition, and so on

[Halil Ibrahim]

first of all thank you guys so much. It was actually a PID Controller System and the cstr is in that for calculating the concertration. I've been working on it for the past 8 hours. I know that seems easy to you but I'm new so, please excuse me. Now I've come up with a solid script I think. It seems the only problem right now is the iteration index in my loop and I couldn't understand it no matter how hard I tried. I've attached the script's final form and the pseudocode for what I have to do plus "Illustrate the system output, reference signal, control signal and evaluation of Kp, Ki , Kd controller parameters for sinusoidal and staircase reference signals as you desired." If you have anything helpful for advancing the script I appreciate that.

ps: I would love to use the built in solver but I have to do it manually with runge kutta 4th order method for my homework but thank you so much for your answer.

[Sam Chak]

Hi @Halil Ibrahim

I'm unfamiliar with your version of the Runge–Kutta 4th-order (RK4) Solver. The RK4 method is typically covered in the Numerical Methods course. In my Numerical Methods course, I solved Mass-Spring-Damper, Pendulum motion, and RLC electrical circuits.

Considering it's merely a homework assignment, I believe this might be too challenging for both of us. Your professor has instructed you to implement a custom RK4 algorithm to solve the CSTR dynamics, (possibly taught at the undergraduate level in Chemical Engineering). Additionally, you are expected to use an Adaptive PID controller, a topic usually encountered at the postgraduate level.

Moreover, you may find it challenging to assess the accuracy of your simulation results obtained with the custom RK4 algorithm unless you perform a comparison and validation. If you opt for the built-in ode45() solver, the code may look something like the following:

%% Simulation settings and call ode45 solver
tspan   = 0:0.01:20;       % simulation time from 0 to 20
x10     = 0.5;
x20     = 0.5;
x30     = 1.5;
x0      = [x10; x20; x30];  % initial values of states x1, x2, x3
[t, x]  = ode45(@cstr, tspan, x0);

%% Plot results
plot(t, x), grid, xlabel Time
legend('x_{1}', 'x_{2}', 'x_{3}')

%% Continuous Stirred-Tank Reactor
function dxdt = cstr(t, x)
    dxdt    = zeros(3, 1);
    
    % Parameters
    d2      = 1;    % this value is hidden somewhere
    u2      = 0;    % maybe the Adaptive PID controller
    
    % Differential equations
    dxdt(1) = 1 - x(1) - 3*x(1) + 0.5*x(2)*x(2);
    dxdt(2) =   - x(2) + 3*x(1) - 0.5*x(2)*x(2) - 1*d2*x(2)*x(2) + u2;
    dxdt(3) =   - x(3)                          + 1*d2*x(2)*x(2);
end

[Halil Ibrahim]

thank you so much for your answer but custom rk4 algorithm is a part of my assignment so its a necessity at this point. Moreover, it doesn't need to be super accurate really.

[Sam Chak]

Hi @Halil Ibrahim,

I'm glad that custom RK4 for Adaptive PID Controller works out for you. However, when I ran the code, it threw the "Index exceeds" error message.

Could you please update the code? This will make your post more meaningful, benefiting users who wish to work with the coding-based Adaptive PID Controller.

clear; clc; close all;

% Initialization of Adaptive PID Control gains
Kp0   = 0.75;
Kp(3) = Kp0;

Ki0   = 0.1;
Ki(3) = Ki0;

Kd0   = 0.25;
Kd(3) = Kd0;

% Control signal saturation parameters
umax  = 1;
umin  = 0;
u(2)  = 0.1;

% CSTR system parameters
Ts    = 0.1;
Da1   = 3;
Da2   = 0.5;
Da3   = 1;

x1(1) = 0.5;
x2(1) = 0.5;
x3(1) = 1.5;

d2    = 1;

for n = 3:100
    % Compute error
    e(n-2) = x1(n-2) - x3(n-2);
    
    % Error-based Adaptive PID Control Signal
    u(n)   = u(n-1) + Kp(n)*(e(n) - e(n-1)) + Ki(n)*e(n) + Kd(n)*(e(n) - 2*e(n-1) + e(n-2));
    
    % Saturate the Control Signal
    if u(n) > umax
        u(n) = umax;
    elseif u(n) < umin
        u(n) = umin;
    end
    
    % Apply Adaptive PID Control Signal to CSTR via custom RK4 algorithm
    k1x1 = Ts*f1_func(x1(n), x2(n));
    k1x2 = Ts*f2_func(x1(n), x2(n), d2, u(n));
    k1x3 = Ts*f3_func(x2(n), x3(n), d2);

    k2x1 = Ts*f1_func(x1(n) + k1x1/2, x2(n) + k1x2/2);
    k2x2 = Ts*f2_func(x1(n) + k1x1/2, x2(n) + k1x2/2, d2, u(n));
    k2x3 = Ts*f3_func(x2(n) + k1x2/2, x3(n) + k1x3/2, d2);

    k3x1 = Ts*f1_func(x1(n) + k2x1/2, x2(n) + k2x2/2);
    k3x2 = Ts*f2_func(x1(n) + k2x1/2, x2(n) + k2x2/2, d2, u(n));
    k3x3 = Ts*f3_func(x2(n) + k2x2/2, x3(n) + k2x3/2, d2);

    k4x1 = Ts*f1_func(x1(n) + k3x1, x2(n) + k3x2);
    k4x2 = Ts*f2_func(x1(n) + k3x1, x2(n) + k3x2, d2, u(n));
    k4x3 = Ts*f3_func(x2(n) + k3x2, x3(n) + k3x3, d2);

    x1(n + 1) = x1(n) + 1/6 * (k1x1 + 2*k2x1 + 2*k3x1 + k4x1);
    x2(n + 1) = x2(n) + 1/6 * (k1x2 + 2*k2x2 + 2*k3x2 + k4x2);
    x3(n + 1) = x3(n) + 1/6 * (k1x3 + 2*k2x3 + 2*k3x3 + k4x3);

    % Update PID gains by adapting to the Error signal
    Kp(n+1) = Kp0 + e(n);
    Ki(n+1) = Ki0 + e(n);
    Kd(n+1) = Kd0 + e(n);

    % Plotting
    figure(1);

    subplot(3,1,1);
    plot(n, x1(n), '-o');
    title('x1');
    hold on;

    subplot(3,1,2);
    plot(n, x2(n), '-o');
    title('x2');
    hold on;

    subplot(3,1,3);
    plot(n, x3(n), '-o');
    title('x3');
    hold on;
end

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

%% Local functions
%  1st CSTR ODE
function dx1 = f1_func(x1, x2)
    t   = 1;
    dx1 = 1 - x1(t) - 3*x1(t) + 0.5*x2(t)*x2(t);
end

%  2nd CSTR ODE
function dx2 = f2_func(x1, x2, d2, u2)
    t   = 1;
    dx2 = -x2(t) + 3*x1(t) - 0.5*x2(t)*x2(t) - 1*d2(t)*x2(t)*x2(t) + u2;
end

%  3rd CSTR ODE
function dx3 = f3_func(x2, x3, d2)
    t   = 1;
    dx3 = -x3(t) + 1*d2(t)*x2(t)*x2(t);
end

[Halil Ibrahim]

Hi @Sam Chak

I'm still trying to fix that index problem you're encounetring when you try to run the code, thats my main problem right now. If I can manage to do that, I will then try to evaluate control parameters using sinusodial and staircase referrance signals, and I'm still open for any help you can give me for that topics. thanks in advance.

[Sam Chak]

Hi @Halil Ibrahim.

I understand that the simulation code for the Adaptive PID controller is not completely fixed. However, I am not familiar with the RK4 version you coded. Now that the issue of passing extra parameters mentioned in your initial question has been resolved, it might be helpful to create a new question specifically addressing the "Index exceeds..." issue.