PI control algorithm manually implemented in matlab

Question

I'm trying to implement a simple script performing a PI control for a cruise control application, but I'm founding some problems with the integral part. Here is my code:

function [] = PI_cruisecontrol()
clc; close all;
t0 = 0; tfinal = 50; dt = 0.001;                % time parameters
r = 10;                                         % reference of 10 m/s
m = 1000;                                       % mass
b = 50;                                         % friction coeff. (depends on v)
yp = zeros(tfinal/dt,1); t = yp;                % initialize speed and time array
Ki = 40;                                        % integrarl constant
Kp = 800;                                       % proportional constant
int = 0;                                        % itinialize int error
% CONTROL LOOP (Forward-Euler integrator is used to solve the ODE)
for i=t0+2:tfinal/dt
    err   = r-yp(i-1);                          % update error
    int   = int+err;                            % integral term
    u     = (Kp*err)+(Ki*int*dt);               % action of control
    yp(i) = yp(i-1)+((-b*yp(i)/m) + (u/m))*dt;  % solve ode for speed  
    t(i)  = t(i)+dt*i;                          % log the time                
end
% Results
figure(1)
plot(t,yp)
title ('Step Response')
xlabel('Time (seconds)')
ylabel('Amplitud')
axis([0 20 0 12])
hold on
reference = ones(tfinal/dt,1)*10;
plot(t,reference,':')
end

And this is how it should be, using predefinided matlab functions:

function [] = PI_cruisecontrol2()
m = 1000;
b = 50;
r = 10;
s = tf('s');
P_cruise = 1/(m*s + b); 
Kp = 800;
Ki = 40;
C = pid(Kp,Ki);
T = feedback(C*P_cruise,1);
t = 0:0.1:20;
step(r*T,t)
axis([0 20 0 12])
end

What am I doing wrong in my code? Thanks!

Answer 1

I managed to fix the problem, working with float variables instead of arrays. Moreover, I added the derivative term (although for this first order problem was not necessary) Here I left the code:

function [] = aFortran_PI()
clc; close all;
r = 10;                         % reference of 10 m/s
m = 1000;                       % mass
b = 50;                         % friction coeff. (depends on v)
yp = 0;                         % init response    
Kp = 800;                       % proportional constant
Ki = 40;                        % proportional constant
Kd = 0;                         % derivative term is not necessary in this problem
previous_error = 0;
integral = 0;
dt = 0.001;
% CONTROL LOOP
for i=1:20000
    error    = r-yp;                        % update error
    integral     = integral + error*dt;     % integral term 
    derivative = (error-previous_error)/dt; % derivative term
    u = Kp*error+Ki*integral+Kd*derivative; % action of control
    yp = yp+(-b*yp/m + u/m)*dt  % solve ode for velocity  

end
end