Index exceeds number of array elements (1)?

    dt = 0.005; N = 200; %(s)
    P3s = 2; P3d = 23.5; %(mmHg)
    V2(1) = .475; %(L)
    V0 = .06; %(L)
    C2 = 0.05; %(L/mmHg)
    R1 = 30; R2 = 30; %(mmHg*s/L)
    P1 = 11; %(mmHg)

    P3(1:120) = P3s;
    P3(121:200) = P3d;
    t = 0:dt:(N-1)*dt; 
    plot(t, P3) 
    xlabel('Time (s)'), ylabel('Pressure (mmHg)')

    for i = 1:N

        P2(i) = (V2(i) - V0) / C2;

        if P2(i) > P3(i); Q3(i) = (P2(i) - P3(i)) / R2; 
            else Q3(i) = 0;
        end

        Q1(i) = (P1 - P2(i)) / R1;

        Q2(i) = Q1(i) - Q3(i);

        P2(i+1) = P2(i) + Q2(i)*dt;
    end

    tp = 0:dt:N*dt; 
    figure
    plot(tp, P2)
    xlabel('Time (s)'), ylabel('P2 (mmHg)')

Error is in line 12 " P2(i) = (V2(i) - V0) / C2;"

What is strange about this problem is that the code runs perfectly on my computer, but when transferred to another and ran, it gives the error message. Any ideas?

[Mark Linne]

I have two nested loops:

for i=1:1:NJ

for f=1:1:NJ

Later in the script I need to do a for loop using:

for L=abs(i-f):1:(i+f)

For obvious reasons, the (i+f) generates the error under discussion. Is there a known way to work around this issue (written by a dinosaur who learned Fortran in 1970, still thinks that way, and would prefer to be able to set array sizes myself at the start of the code)?

Answer 1

[Aquatris]

You only define V2(1) in your code and do not increase the size of the V2 variable. So I think you are missing a code in your for loop which assigns values to V2(i+1).

It probably runs in your computer because your workspace includes a V2 variable with correct size from another script so the code does not yell at you.

[Riley Vargas]

Thank you! It worked, I just had to add a line of code that defined V2(i+1). You were right.

[Lykha ahmed]

could you tell me how you added?

[parag gupta]

how did u make the change?

[Aquatris]

In his code he never defined V2(2) so he added a line in his for loop, something like;

for i = 1:N
      ...
      V2(i+1) = *equation*
end

Answer 2

[sura naji]

clc
clear
close all
% Problem Statement
Npar = 17;
VarLow=0.1;
VarHigh = 35;
%BBBC parameters
N=50;       %number of candidates
MaxIter=100;  %number of iterations
% initialize a random value as best value
XBest = rand(1,Npar).* (VarHigh - VarLow) + VarLow;
FBest=fitnessFunc(XBest);
GB=FBest;
t = cputime;
%intialize solutions and memory
X = zeros(N, Npar);
F = zeros(N, 1);
for ii = 1:N
    X(ii,:) = rand(1,Npar).* (VarHigh - VarLow) + VarLow;
    
    % calculate the fitness of solutions
    F(ii) = fitnessFunc(X(ii,:));
end
%Main Loop
for it=1:MaxIter
    %Find the centre of mass 
    %-----------------------
    %numerator term
    num=zeros(1,Npar);
    for ii=1:N
        for jj=1:Npar
            num(jj)=num(jj)+(X(ii,jj)/F(ii));
        end
    end
    %denominator term
    den=sum(1./F);
    %centre of mass
    Xc=num/den; 
    %generate new solutions
    %----------------------
    for ii=1:N
        %new solution from centre of mass
        for jj=2:Npar      
            New=X(ii,:);
            New(jj)=Xc(jj)+((VarHigh(jj)*rand)/it^2);
            
        end
        %boundary constraints
        New=limiter(New,VarHigh,VarLow);
        %new fitness
        newFit=fitnessFunc(New);
        %check whether the solution is better than previous solution
        if newFit<F(ii)
            X(ii,:)=New;
            F(ii)=newFit;
            if F(ii)<FBest
                XBest=X(ii,:);
                FBest=F(ii);   
            end
        end
    end
    % store the best value in each iteration
    GB=[GB FBest];
end
%Main Loop
for it=1:MaxIter
    %Find the centre of mass 
    %-----------------------
    %numerator term
    num=zeros(1,Npar);
    for ii=1:N
        for jj=1:Npar
            num(jj)=num(jj)+(X(ii,jj)/F(ii));
        end
    end
    %denominator term
    den=sum(1./F);
    %centre of mass
    Xc=num/den; 
    %generate new solutions
    %----------------------
    for ii=1:N
        %new solution from centre of mass
        for jj=2:Npar      
            New=X(ii,:);
            New(jj)=Xc(jj)+((VarHigh(jj)*rand)/it^2);
            
        end

hi please can you help me because l have the same problem and the message related to [ New(jj)=Xc(jj)+((VarHigh(jj)*rand)/it^2);] [index exceeds number of array element(1)]?

[Aquatris]

VarHigh variable is a scalar not vector since you only defined it with "VarHigh = 35;" line.

And in the part that gives you the error, you are trying to reach "VarHigh(2)" which does not exist.

Answer 3

[BAMIDELE JEGEDE]

I have the issue with my code

clear, clc
x = 0:0.01:pi
y = myfunc(x)
plot(x,y)

avg_y = y(1:length(x)-1) + diff(y)/2;
A = sum(diff(x).*avg_yin )

at 
avg_y = y(1:length(x)-1) + diff(y)/2;

I honestly don't know what I am doing wrong and it's becoming frustrating

[Aquatris]

I replaced myfunc with sin function and it does not give any error at line "avg_y = y(1:length(x)-1) + diff(y)/2;"

However "A = sum(diff(x).*avg_yin )" this one gives error since avg_yin is not defined in the code you gave.

What is the exact error message you are getting?

[Farshid R]

Hi

Good time

I wrote this code but it gives an error

Please help me

thank you

n=100;
u1=[0,0]';
X1=[-4,-2,0]';



            % p=data.PredictionHorizon;
        a=0.9;h=0.9;
        cp1=1;cp2=1;cp3=1;

            % c1=0;c2=0;c3=0;
            for j=1:n
                c1(j)=(1-(1+a)/j)*cp1;
                c2(j)=(1-(1+a)/j)*cp2;
                c3(j)=(1-(1+a)/j)*cp3;
                cp1=c1(j); cp2=c2(j); cp3=c3(j);
            end
            % initial conditions setting:
        v1(1)=u1(1);
        w1(1)=u1(2);
        x1(1)=X1(1); y1(1)=X1(2); z1(1)=X1(3);
            % calculation of phase portraits /numerical solution/:
            for i=2:n
                x1(i)=h*cos(z1(i-1))*v1(i-1) - memo(x1, c1, i);
                y1(i)=h*sin(z1(i-1))*v1(i-1)-memo(y1, c2, i);
                z1(i)=h*w1(i-1)-memo(z1, c3, i) ;
            end
%%
function [yo] = memo(r, c, k)
%
temp = 0;
for j=1:k-1
   temp = temp + c(j)*r(k-j);
end
yo = temp;
%
%%%%% error

Index exceeds the number of array elements (1).

Error in exocstrstateFcnCT1 (line 28)

x1(i)=h*cos(z1(i-1))*v1(i-1) - memo(x1, c1, i);

[Walter Roberson]

        v1(1)=u1(1);

You initialize one v1 value

            for i=2:n
                x1(i)=h*cos(z1(i-1))*v1(i-1) - memo(x1, c1, i);

When i = 3, you access v1(3-1) -> v1(2) . But that does not exist.

                y1(i)=h*sin(z1(i-1))*v1(i-1)-memo(y1, c2, i);
                z1(i)=h*w1(i-1)-memo(z1, c3, i) ;
            end

You keep extending x1 and y1 and z1 in that loop, but you do not extend v1 or w1

Answer 4

[Sara Babaie]

I need some help with my code...I am getting the same error:

format long
ep=0.8; %Radiative emissivity of the pipe surface
stfb=5.67e-08; %Stefan-Boltzmann constant
q=500; %heat rate
tair=270; %air temperature
tsurr=270; %surrounding temperature
h=10; %convective heat coefficient
d=0.1; %diameter of the pipe
L=0.5; %length of pipe
p=pi;
ft=@(t84)q-(p*d*L(h*(t84-tair)+(ep*stfb(t84^4-tsurr^4))))
a=ft(334)

[Walter Roberson]

stfb=5.67e-08; %Stefan-Boltzmann constant

That is a scalar.

stfb(t84^4-tsurr^4)

and that tries to index the scalar at a large integer.

MATLAB has absolutely NO implied multiplication. Not anywhere. Not even in the internal workings of the symbolic computing engine that few dare look at.

[Farshid R]

This error is related to the dimensions and the number of elements is more.

[Walter Roberson]

In this case, Sara's error is in forgetting a multiplication symbol. stfb*(t84^4-tsurr^4)

Answer 5

[Simon Hudon-Labonté]

function [xsol]=NewtonRaphson(Fon,DerFon,SolEst,Erreur,jmax)
iflag=0;
for j=1:jmax
xsoli=SolEst-(Fon(SolEst)/DerFon(SolEst));
if abs((xsoli-SolEst)/SolEst)<Erreur
xsol=xsoli;
iflag=1;
break
end
SolEst=xsoli;
end
if j==jmax && iflag~=1
fprintf('Aucune solution après %i itérations.\n',jmax)
xsol=('Aucune réponse');
end

Hello i am getting the same error as the others in the comment but even with the different answers I cannot find my error. Thank you!

Simon

[Farshid R]

This error is related to the dimensions and the number of elements is more.

[Walter Roberson]

You would get that error if the first parameter was an array instead of a function handle

[Binyam Sime]

Hello!!

I need some little correction (If any) for the code blow!

If any one help me, i have a greate gratittude!!!!!

A = imread('C:\Users\pc_2\Desktop\bs\Im_one.jpg');
r=size(A,1);
c=size(A,2);
A=double(A);
ah= uint8(zeros(r,c));
n=r*c;
f=zeros(256,1);
pdf=zeros(256,1);
cdf=zeros(256,1);
cum=zeros(256,1);
out=zeros(256,1);
for i=1:r
    for j=1:c
        value =A(i,j);
        f(value+1)=f(value+1)+1;
        pdf(value+1) =f(value+1)/n;
    end
end
sum=0;L=255;

for i=1:size(pdf)
    sum =sum +freq(i);
    cum(i)=sum;
    cdf(i)=cum(i)/n;
    out(i) =round(cdf(i)*L);
end
for i=1:r
    for j=1:c
        ah(i,j)= out(A(i,j)+1);
    end
end
figure,imshow(ah);
he= histeq(A);
figure,imshow(he);

when run it says as follws:-

Index exceeds the number of array elements (1).

Thank you for your help!

[Walter Roberson]

    sum =sum +freq(i);

You do not define freq() in the code.

Answer 6

[Wu Changjin Wu]

can some one help me, I got the same error !

% Read in image.
grayImage = imread('FFT.JPG');
[rows, columns, numberOfColorChannels] = size(grayImage)
if numberOfColorChannels > 1
  grayImage = rgb2gray(grayImage);
end
% Display original grayscale image.
subplot(2, 3, 1);
imshow(grayImage)
axis on;
title('Original Gray Scale Image', 'FontSize', fontSize)
% Enlarge figure to full screen.
set(gcf, 'units','normalized','outerposition',[0 0 1 1]);
% Perform 2D FFTs
fftOriginal = fft2(double(grayImage));
% Move center from (1,1) to (129, 129) (the middle of the matrix).
shiftedFFT = fftshift(fftOriginal);
subplot(2, 3, 2);
scaledFFTr = 255 * mat2gray(real(shiftedFFT));
imshow(log(scaledFFTr), []);
title('Log of Real Part of Spectrum', 'FontSize', fontSize)
subplot(2, 3, 3);
scaledFFTi = mat2gray(imag(shiftedFFT));
imshow(log(scaledFFTi), []);
axis on;
title('Log of Imaginary Part of Spectrum', 'FontSize', fontSize)
% Display magnitude and phase of 2D FFTs
subplot(2, 3, 4);
shiftedFFTMagnitude = abs(shiftedFFT);
imshow(log(abs(shiftedFFTMagnitude)),[]);
axis on;
colormap gray
title('Log Magnitude of Spectrum', 'FontSize', fontSize)
% Get average
midRow = rows/2+1
midCol = columns/2+1
maxRadius = ceil(sqrt(129^2 + 129^2))
radialProfile = zeros(maxRadius, 1);
count = zeros(maxRadius, 1);
for col = 1 : columns
  for row = 1 : rows
    radius = sqrt((row - midRow) ^ 2 + (col - midCol) ^ 2);
    thisIndex = ceil(radius) + 1;
    radialProfile(thisIndex) = radialProfile(thisIndex) + shiftedFFTMagnitude(row, col);
    count(thisIndex) = count(thisIndex) + 1;
  end
end
radialProfile = radialProfile./ count;
subplot(2, 3, 5:6);
plot(radialProfile, 'b-', 'LineWidth', 2);
grid on;
title('Average Radial Profile of Spectrum', 'FontSize', fontSize)

[Walter Roberson]

fontSize = 12;

% Read in image.
filename = 'FFT.JPG';
if isunix()
  filename = 'flamingos.jpg';
end

grayImage = imread(filename);
[rows, columns, numberOfColorChannels] = size(grayImage);
if numberOfColorChannels > 1
  grayImage = rgb2gray(grayImage);
end
% Display original grayscale image.
subplot(2, 3, 1);
imshow(grayImage)
axis on;
title('Original Gray Scale Image', 'FontSize', fontSize)
% Enlarge figure to full screen.
set(gcf, 'units','normalized','outerposition',[0 0 1 1]);
% Perform 2D FFTs
fftOriginal = fft2(double(grayImage));
% Move center from (1,1) to (129, 129) (the middle of the matrix).
shiftedFFT = fftshift(fftOriginal);
subplot(2, 3, 2);
scaledFFTr = 255 * mat2gray(real(shiftedFFT));
imshow(log(scaledFFTr), []);
title('Log of Real Part of Spectrum', 'FontSize', fontSize)
subplot(2, 3, 3);
scaledFFTi = mat2gray(imag(shiftedFFT));
imshow(log(scaledFFTi), []);
axis on;
title('Log of Imaginary Part of Spectrum', 'FontSize', fontSize)
% Display magnitude and phase of 2D FFTs
subplot(2, 3, 4);
shiftedFFTMagnitude = abs(shiftedFFT);
imshow(log(abs(shiftedFFTMagnitude)),[]);
axis on;
colormap gray
title('Log Magnitude of Spectrum', 'FontSize', fontSize)
% Get average
midRow = rows/2+1;
midCol = columns/2+1;
%maxRadius = ceil(sqrt(129^2 + 129^2))
maxRadius = ceil(sqrt(midRow.^2 + midCol.^2));

radialProfile = zeros(maxRadius, 1);
count = zeros(maxRadius, 1);
for col = 1 : columns
  for row = 1 : rows
    radius = sqrt((row - midRow) ^ 2 + (col - midCol) ^ 2);
    thisIndex = ceil(radius) + 1;
    radialProfile(thisIndex) = radialProfile(thisIndex) + shiftedFFTMagnitude(row, col);
    count(thisIndex) = count(thisIndex) + 1;
  end
end
radialProfile = radialProfile./ count;
subplot(2, 3, 5:6);
plot(radialProfile, 'b-', 'LineWidth', 2);
grid on;
title('Average Radial Profile of Spectrum', 'FontSize', fontSize)

Answer 7

[ELIZABETH ESPARZA]

Hi, i have the same error, i need help :(

clear all 

%incremento en x
h=0.1;

%funcion 1
f = @(t,a,b,u,v,w) -5*a*u+2*a-4*b;

%funcion 2
g = @(t,a,b,u,v,w) -a*b*v-a+6*b;

%funcion 3
d = @(t,a,b,u,v,w) b*w-a;


%Condiciones Iniciales
ti = 0;

ai = 0; 

bi = 0;

ui = 0;

vi = 0;

wi = 00;

%Iteraciones
n = 5;

%Función
[t,a,b,u,v,w] =RungeFG3(f,g,d,n,h,ti,ai,bi,ui,vi,wi)

%Graficación
plot3(t,a,b);grid on;
legend('Runge K4')
title('Gráfica')
xlabel('t')
ylabel('a')
zlabel('b')

tabla1= table (t',a',b');
tabla1.Properties.VariableNames = {'t','a','b'}

function [t,a,b,u,v,w] =RungeFG3(f,g,d,n,h,t0,a0,b0,u0,v0,w0)

t(1) = t0;
a(1) = a0;
b(1) = b0;
u(1) = u0;
v(1) = v0;
w(1) = w0;

for i=1:1:n
t(i+1) = t(i) + h;

k1 = h*f(t(i),a(i),b(i),u(i),v(i),w(i));

l1 = h*g(t(i),a(i),b(i),u(i),v(i),w(i));
e1 = h*d(t(i),a(i),b(i),u(i),v(i),w(i));

m1 = h*u(i);
n1 = h*v(i);
p1 = h*w(i);

k2 = h*f(t(i)+h/2,a(i)+m1/2,b(i)+n1/2,u(1)+k1/2,v(i)+l1/2,w(i)+l1/2);
l2 = h*g(t(i)+h/2,a(i)+m1/2,b(i)+n1/2,u(1)+k1/2,v(i)+l1/2,w(i)+l1/2);
e2 = h*d(t(i)+h/2,a(i)+m1/2,b(i)m2 = h*(u(i)+k1/2);
n2 = h*(v(i)+l1/2);
p2 = h*(w(i)+p1/2);

k3 = h*f(t(i)+h/2,a(i)+m2/2,b(i)+n2/2,u(1)+k2/2,v(i)+l2/2,w(i)+l2/2);
l3 = h*g(t(i)+h/2,a(i)+m2/2,b(i)+n2/2,u(1)+k2/2,v(i)+l2/2,w(i)+l2/2);
e3 = h*d(t(i)+h/2,a(i)+m2/2,b(i)+n2/2,u(1)+k2/2,v(i)+l2/2,w(i)+l2/2);

m3 = h*(u(i)+k2/2);
n3 = h*(v(i)+l2/2);
p3 = h*(w(i)+p2/2);

k4 = h*f(t(i)+h,a(i)+m3,b(i)+n3,u(1)+k3,v(i)+l3,w(i)+l3);
l4 = h*g(t(i)+h,a(i)+m3,b(i)+n3,u(1)+k3,v(i)+l3,w(i)+l3);
e4 = h*d(t(i)+h,a(i)+m3,b(i)+n3,u(1)+k3,v(i)+l3,w(i)+l3);

m4 = h*(u(i)+k3);
n4 = h*(v(i)+l3);
p4 = h*(w(i)+p3);


u(i+1) = u(i)+(k1+2*k2+2*k3+k4)/6;
v(i+1) = v(i)+(l1+2*l2+2*l3+l4)/6;
p(i+1) = w(i)+(p1+2*p2+2*p3+p4)/6;+n1/2,u(1)+k1/2,v(i)+l1/2,w(i)+l1/2);

a(i+1) = a(i)+(m1+2*m2+2*m3+m4)/6;
b(i+1) = b(i)+(n1+2*n2+2*n3+n4)/6;

end
end 

[Walter Roberson]

You have a different error in the code you posted

e2 = h*d(t(i)+h/2,a(i)+m1/2,b(i)m2 = h*(u(i)+k1/2);
        1 2 1      2 1       2 1   ?

You do not have an operator between the b(i) and the m2, and after the m2 you have an = but you are already within a = .

I suspect there should be a ); after the b(i)

[Walter Roberson]

p(i+1) = w(i)+(p1+2*p2+2*p3+p4)/6;+n1/2,u(1)+k1/2,v(i)+l1/2,w(i)+l1/2);

The semi-colon marks the end of the assigment statement.

After that you have +n1/2 as a statement. That is unusual as a statement, but legal.

After that you have u(1)+k1/2 as a statement. That is legal.

After that you have v(1)+l1/2 as a statement. That is legal.

After that you have w(i)+l1/2) as a statement. That is not legal as it has an extra ) that does not match any (

It looks to me as if perhaps everything after the semi-colon was intended to be part of the assignment to e2.

[ELIZABETH ESPARZA]

[ELIZABETH ESPARZA]

oh well, but the mistake is in line 58

[Aquatris]

That index problem happens because of variable w. You never assigned values to w(i+1) in your code.

Answer 8

[emre bahçeci]

Hello,i am trying to face same error.Here is my code:

T=1.5
t=0:0.001:T
m=1;b=10;k=100;
K=[1 0 0 0 0 0;
   0 1 0 0 0 0;
   0 0 2 0 0 0;
   1 T T^2 T^3 T^4 T^5;
   0 1 2*T 3*T^2 4*T^3 5*T^4;
   0 0 2 6*T 12*T^2 20*T^3]
A=inv(K)*[pi/3;0;0;5*pi/3;0;0]
syms theta(x) 
theta_dot=diff(theta,x)
theta_dot_dot=diff(theta_dot,x)
Q=m*theta_dot_dot+b*theta_dot+k*theta
for i=1:length(t)
    teta(i)=A(1)+A(2)*t(i)+A(3)*t(i)^2+A(4)*t(i)^3+A(5)*t(i)^4+A(6)*t(i)^5
    tetah(i)=A(2)+2*A(3)*t(i)+3*A(4)*t(i)^2+4*A(5)*t(i)^3+5*A(6)*t(i)^4
    tetai(i)=2*A(3)+6*A(4)*t(i)+12*A(5)*t(i)^2+20*A(6)*t(i)^3
    Q_new(i)=subs(Q,[theta_dot_dot theta_dot theta],[tetai tetah teta])
end

After that i have encountered same error:

Index exceeds the number of array elements (1).

Error in sym/subs>normalize (line 212)

Y{i} = fun2sym(Y{i},argnames(X{i}));

Error in sym/subs>mupadsubs (line 166)

[X2,Y2,symX,symY] = normalize(X,Y); %#ok

Error in sym/subs (line 154)

G = mupadsubs(F,X,Y);

Error in symbolic_in_for_loop (line 19)

Q_new(i)=subs(Q,[theta_dot_dot theta_dot theta],[tetai tetah teta])

Can anybody help me?

Thanks

[Walter Roberson]

I had to make T smaller for demonstration purposes

format long g
T=.015;
t=0:0.001:T;
m=1;b=10;k=100;
K=[1 0 0 0 0 0;
   0 1 0 0 0 0;
   0 0 2 0 0 0;
   1 T T^2 T^3 T^4 T^5;
   0 1 2*T 3*T^2 4*T^3 5*T^4;
   0 0 2 6*T 12*T^2 20*T^3]

K = 6×6

1 0 0 0 0 0 0 1 0 0 0 0 0 0 2 0 0 0 1 0.015 0.000225 3.375e-06 5.0625e-08 7.59375e-10 0 1 0.03 0.000675 1.35e-05 2.53125e-07 0 0 2 0.09 0.0027 6.75e-05

A=inv(K)*[pi/3;0;0;5*pi/3;0;0]

A = 6×1

1.0471975511966 0 0 12411230.2364042 -1241123023.64042 33096613963.7446

syms theta(x) 
theta_dot=diff(theta,x)

theta_dot(x) = 

theta_dot_dot=diff(theta_dot,x)

theta_dot_dot(x) = 

Q=m*theta_dot_dot+b*theta_dot+k*theta

Q(x) = 

for i=1:length(t)
    teta(i)=A(1)+A(2)*t(i)+A(3)*t(i)^2+A(4)*t(i)^3+A(5)*t(i)^4+A(6)*t(i)^5;
    tetah(i)=A(2)+2*A(3)*t(i)+3*A(4)*t(i)^2+4*A(5)*t(i)^3+5*A(6)*t(i)^4;
    tetai(i)=2*A(3)+6*A(4)*t(i)+12*A(5)*t(i)^2+20*A(6)*t(i)^3;
    Q_new(i,1)=subs(Q(x),{theta_dot_dot theta_dot theta},{tetai(i) tetah(i) teta(i)});
end

Q_new

Q_new = 

[emre bahçeci]

Thanks a lot

Answer 9

[Anastasiya Moiseeva]

Please, can someone help me?

I got the same error in my m-file

Error in file (line 16)

yp = simout(b)

% Data

u1 = 40;

u2 = 60;

u = 50;

y1 = 0;

y2 = 60;

y3 = 40;

mui1 = 0;

mui2 = 100;

time1 = 0;

%Calculate parameters

% inflection point

[a, b] = max(yout)

tp = tout(b)

% simout = round(simout) % rounding function

yp = simout(b)

% Time time2

d = find(simout==60)

time2 = tout(d)

% Time time3

e = find(simout==40)

time3 = tout(e)

%Square Sy2

u = -yout(e)

Sy2 = ((u1-y1)*(u1-y1))/u

%Square Sy1

Sy1 = ((sum(simout))*1.65) %/ it is necessary to change the scale in simout

Smui = (mui2-mui1)*(time2-time1)

%Find Коб

Kob = (Sy1+Sy2)/Smui

%Find b

b = (yp-y1)/(Kob*(mui2-mui1))

% Find x

if (b < 0.14)

c1 = 0.2604;

c2 = -0.07986;

c3 = -0.0203;

x = c2+c3/(b-c1)

elseif (0.14 <= b) && (b < 0.26)

c1 = 0.2993;

c2 = -0.1076;

c3 = -0.03128;

x = c2+c3/(b-c1)

else

x = 0.6881

end

% Find z

z = x^(x/(1-x))

% Find Tob

Tob = (Kob*(mui2-mui1))/a

% Find the parameters of the object model (time constants ... object models

Tmob2 = z*Tob

Tmob1 = x*Tmob2

taumob = tp-Tmob2*[-log(z)]

%Object Model Parameters

taumob = taumob;

Tmob1 = Tmob1;

Tmob2 = Tmob2;

Kmob = Kob;

M = 1.4;

%

%Values .......................Rrs.op....Frs.op.........

Rs.op = 1.4; %sqrt(b^2 + c^2) 1.3

Gs.op = -90; %atand(-c/b)

Fs.op = Gs.op*pi/180

Ws.op = Rs.op*exp(j*Fs.op)

%====-----------------------

Wrs.op = Ws.op/(1-Ws.op)

%====Wrs.o

%aa = real(Wrs.op)

%bb = imag(Wrs.op)

Rrs.op = abs(Wrs.op)

Frs.op = angle(Wrs.op)

%---------------------------------------

%Determination of the optimal values of the aprha parameters *****

n = Tmob2/Tmob1

beta = taumob/Tmob1

a01 = -4;

b01 = 0.38;

c01 = 0.2;

a16 = -3.714286;

b16 = 0.1817143;

c16 = 0.5576327;

az4 = -0.2613139;

bz4 = 0.2176277;

cz4 = 0.0677002;

y01 = b01 + (c01/(n-a01));

y16 = b16 + (c16/(n-a16));

z4 = bz4 + (cz4/(beta-az4));

bbb = -6.622517;

ccc = -2.5821192;

d4 = bbb*z4 -ccc;

d4 = bbb*z4 - ccc;

anpha = y01*(1-d4)+y16*d4

%Determination of optimal parameter values ТТi.op ****

a011 = -0.10738;

b011 = 1.25235;

c011 = 0.60646;

a161 = -2.20225;

b161 = 1.60899;

c161 = 8.93751;

az41 = -1.4;

bz41 = 4.7;

cz41 = -4.95;

bbb1 = 0.60606;

ccc1 = 0.84848;

y011 = b011 + (c011/(n-a011));

y161 = b161 + (c161/(n-a161));

z41 = bz41 + (cz41/(beta - az41));

d41 = bbb1*z41 - ccc1;

TTi.op = y011*(1-d41) + y161*d41

%----------------------------------------

Kf = 8;

a_Ti = 2*pi/TTi.op;

a_Td = a_Ti*anpha;

a_Tf = a_Td/Kf;

Wf_op = 1/((j*a_Tf + 1)*(j*a_Tf + 1));

Rf_op = abs(Wf_op)

Ff_op = angle(Wf_op)

Wr_op = 1 + (1/(j*a_Ti)) + j*a_Td*Wf_op;

Rr.op = abs(Wr_op)

Fr.op = angle(Wr_op)

Fob.op = Frs.op

% The optimal value of the dimensionless frequency is determined

x1 = 1.5;

Gx = beta*x1 + atan(x1) + atan(x1*n) + Fob.op;

GGx = beta + 1/(x1^2 +1) + n/((n^2)*(x1^2) +1);

finish = abs(Gx/GGx);

for i = 1:100

if (finish >= 0.01)

x1 = x1 - Gx/GGx;

Gx = beta*x1 + atan(x1) + atan(x1*n) + Fob.op;

GGx = beta + 1/(x1^2 +1) + n/((n^2)*(x1^2) +1);

finish = abs(Gx/GGx);

else

x0 = x1;

Om_op = x1

break

end

i = i+1;

end

Ks.op = (Rrs.op*sqrt((Om_op^2 + 1)*((Om_op^2)*n^2 + 1)))/Rr.op

Kr.op = Ks.op/Kmob

T0 = (Tmob1*2*pi)/Om_op;

Ti.op = T0/TTi.op

Td = Ti.op*anpha

%-Controller Parameters

Kp = Kr.op

Ki = Kp/Ti.op

Kd = Kp*Td

Tf = Td/Kf

%----------------------------------------

[Walter Roberson]

simout is not defined in your code. If you assigned

simout = sim(SOME Simulink stuff)

then you probably configured Simulink to return "unified" output, which is a scalar struct with fields named after the variables being returned.

Answer 10

[Ivan Bullier Ananfack Tsakem]

Index exceeds the number of

array elements (1).

Error in Trajectory/show (line

73)

carPositionOnCenterline =

obj.nearest_points(car_states(1),

car_states(2));

Error in Main (line 78)

myTrajectory.show(Lambo,

file_path)

>>

[Ivan Bullier Ananfack Tsakem]

classdef Trajectory < handle

properties

WP_World_Coordinate;

nump ;

Nearst_Points;

Nearst_Points_C;

cte;

show_nearest = true;

show_fit = true;

Coeff_fit;

centerline;

end

methods

%% Constructor

function obj = Trajectory(way_points, nump)

% Setze obj.nump auf einen sinnvollen Standardwert, wenn nump nicht angegeben ist

if nargin < 2 || isempty(nump) || nump <= 0

obj.nump = 400; % oder einen anderen geeigneten Standardwert

else

obj.nump = nump;

end

obj.WP_World_Coordinate = way_points;

obj.Nearst_Points = zeros(obj.nump, 2);

obj.Nearst_Points_C = zeros(obj.nump, 2);

obj.cte = 0;

obj.Coeff_fit = [0 0 0];

obj.centerline = []; % Hinzugefügt: Initialisiere Centerline als leere Matrix

end

%% Show waypoints

function show(obj, Car, file_path)

% Read Pylon data from file

%pylons_data = readmatrix(file_path);

pylons_data = readtable(file_path, 'Delimiter', ',', 'Format', '%s%f%f');

disp(pylons_data);

% Überprüfen Sie die Spaltennamen

disp(class(pylons_data));

%disp(pylons_data.Properties.VariableNames);

% Überprüfen Sie die ersten paar Zeilen der Daten

blue_pylons = pylons_data(strcmp(pylons_data.Var1, 'blue'), :);

yellow_pylons = pylons_data(strcmp(pylons_data.Var2, 'yellow'), :);

if height(blue_pylons) < obj.nump || height(yellow_pylons) < obj.nump

disp('Warnung: Die Anzahl der Pylone ist nicht ausreichend.');

obj.nump = max(obj.nump, min(height(blue_pylons), height(yellow_pylons)));

disp(['Angepasste Anzahl von obj.nump: ', num2str(obj.nump)]);

end

blue_pylons_length = min(height(blue_pylons), obj.nump);

yellow_pylons_length = min(height(yellow_pylons), obj.nump);

figure; % Neue Figur erstellen

scatter(blue_pylons.Var2, blue_pylons.Var3, 'bo', 'LineWidth', 1.5);

hold on

scatter(yellow_pylons.Var2, yellow_pylons.Var3, 'yo', 'LineWidth', 1.5);

hold on

% Überprüfen Sie, ob die Centerline leer oder gültig ist

if ~isempty(obj.centerline) && size(obj.centerline, 1) > 1

% Plot Centerline

obj.draw_centerline();

else

% Die Centerline ist leer oder ungültig, geben Sie eine Warnung aus

warning('Die Centerline ist leer oder nicht korrekt zugewiesen.');

end

% Hier den Zustand des Autos abrufen und die Variable carPositionOnCenterline aktualisieren

car_states = Car.state_unpack();

carPositionOnCenterline = obj.nearest_points(car_states(1), car_states(2));

if ~isempty(carPositionOnCenterline)

car_x_pid = carPositionOnCenterline(1);

car_y_pid = carPositionOnCenterline(2);

else

error('nearest_points returns an empty array.');

end

if numel(carPositionOnCenterline) >= 2

car_x_pid = carPositionOnCenterline(1);

car_y_pid = carPositionOnCenterline(2);

else

disp('Warnung: Ungültige Größe von carPositionOnCenterline.');

car_x_pid = 0;

car_y_pid = 0;

end

% Restlicher Code bleibt unverändert

disp(obj.centerline);

disp(car_x_pid);

disp(car_y_pid);

dist = sqrt((obj.centerline(:, 1) - car_x_pid).^2 + (obj.centerline(:, 2) - car_y_pid).^2);

car.states(1:2) = [car_x_pid, car_y_pid];

car.show();

hold off;

end

function centerline = find_center_points(obj, file_path)

% Read pylon data from file

pylons_data = readtable(file_path, 'Delimiter', ',', 'Format', '%s%f%f');

% Überprüfen, ob pylons_data nicht leer ist

if isempty(pylons_data)

warning('Pylon data is empty. Check the data file.');

centerline = []; % oder einen Standardwert, je nach Bedarf

return;

end

% Verwende die Methode aus dem Kommentar für die Berechnung der Centerline

% Verwende curly brace indexing, um direkt auf die Daten zuzugreifen

pylon_positions = pylons_data{:, 2:3};

centerline = (pylon_positions(1:end-1, :) + pylon_positions(2:end, :)) / 2;

% Speichere die berechnete Centerline in der Eigenschaft der Klasse

obj.centerline = centerline;

end

function draw_centerline(obj)

% Plot Centerline

plot(obj.centerline(:, 1), obj.centerline(:, 2), 'r-', 'LineWidth', 2);

num_points = size(obj.centerline, 1);

% Füge eine Legende für die Centerline hinzu

legend(['Number of Points: ', num2str(num_points)], 'Location', 'best');

end

%% Transform the point to car coordinate system

function W_to_Car_Coordinate_system(obj, Car)

[x, y, phi, ~] = Car.state_unpack;

% Korrigiere den Tippfehler im Variablennamen

translate = obj.Nearst_Points - [x, y];

obj.Nearst_Points_C = ([cos(phi) sin(phi); -sin(phi) cos(phi)] * translate')';

end

%% This function finds the nearest points to the car

function nearest_center = nearest_points(obj, car_x_pid, car_y_pid, centerline)

% Überprüfen, ob centerline nicht leer ist

if isempty(obj.centerline)

warning('The Centerline is empty. Check the pylon data.');

nearest_center = []; % oder einen Standardwert, je nach Bedarf

return;

end

% Alternative to knnsearch without Statistics and Machine Learning Toolbox

centerline = obj.centerline;

distances = sqrt((obj.centerline(:, 1) - car_x_pid).^2 + (obj.centerline(:, 2) - car_y_pid).^2);

[~, indices] = sort(distances);

nearest_indices = knnsearch(obj.centerline, [car_x_pid, car_y_pid], 'K', obj.nump);

nearest_center = obj.centerline(nearest_indices, :);

end

function poly_fit(obj, Car)

obj.W_to_Car_Coordinate_system(Car);

x = obj.Nearst_Points_C(:, 1);

y = obj.Nearst_Points_C(:, 2);

% Zentrieren und skalieren der Daten

x_centered = x - mean(x);

x_scaled = x_centered / std(x_centered);

p = polyfit(x_scaled, y, 2);

obj.Coeff_fit = [p(1) p(2) p(3)];

end

%% This function computes the track error which will be used later by the PID controller

function compute_error(obj, car_x_pid, car_y_pid)

nearest_center = obj.nearest_points(car_x_pid, car_y_pid);

obj.cte = norm([car_x_pid, car_y_pid] - nearest_center);

end

function centerline = getCenterline(obj)

centerline = obj.centerline;

end

end

end

%% Transform the point to car coordinate system

% function W_points = Car_to_W_Coordinate_system(C_Points, Car)

% [x, y, phi, ~] = Car.state_unpack;

% Rotate = ([cos(phi) -sin(phi); sin(phi) cos(phi)] * C_Points')';

% W_points = Rotate + [x, y];

% end

%

[Ivan Bullier Ananfack Tsakem]

close all

clc

% Plot options

x0 = 0;

y0 = 0;

gold = (1 + sqrt(5)) / 2;

width = 10;

height = 4;

FontSize = 14;

FontSize_axis = 16;

font_weight = 'bold';

% Create your object

initial_states = [0 -225 0 0];

initial_inputs = [0 12];

Lambo = Car(initial_states, initial_inputs);

N = 200; % Number of time steps for the simulation.

win = 120;

way_points = generate_trajectory(0.05);

file_path = 'C:\Users\ivant\OneDrive\Bureau\Neuer Ordner\Self_Driving_Car_Simulation_In_MATLAB-master\Self_Driving_Car_Simulation_In_MATLAB-master\spline_dataset_400p.csv';

% Create a Trajectory object

myTrajectory = Trajectory(way_points);

% Find Centerline

centerline = myTrajectory.find_center_points(file_path);

% Check if the Centerline is not empty

if isempty(centerline)

% The Centerline is empty, issue a warning

warning('The Centerline is empty. Check the pylon data.');

return; % Exit the script if the Centerline is empty

end

% Gif file configuration

filename = 'trajectory_simulation.gif';

delay_time = 0.7;

h = figure;

% Preallocate arrays for speed

x = zeros(N, 1);

y = zeros(N, 1);

theta = zeros(N, 1);

v = zeros(N, 1);

% Vectorize the loop to avoid repeated function calls

[x, y, theta, v] = Lambo.state_unpack;

x(isnan(x)) = 0;

y(isnan(y)) = 0;

% Find the nearest points on the Centerline to the car

%nearest_center = myTrajectory.nearest_points(x, y, centerline);

nearest_center = myTrajectory.nearest_points(Lambo.states(1), Lambo.states(2));

% Fit a polynomial to the nearest points

myTrajectory.poly_fit(Lambo);

% Compute the cross-track error

myTrajectory.compute_error(Lambo.states(1), Lambo.states(2));

% Apply the PID controller

Lambo.PID_Controller(myTrajectory.cte);

% Control the inputs

Lambo.control_inputs(1);

% Update the state

Lambo.update_state;

% Plot the car, the Centerline, and the nearest points

%Lambo.show;

myTrajectory.draw_centerline();

myTrajectory.show(Lambo, file_path)

carPositionOncenterline = myTrajectory.nearest_points(Lambo.getStates(1), Lambo.getStates(2));

carPositionOncenterline = obj.positionOnCenterline(Car);

myTrajectory.centerline = myTrajectory.centerline(~any(isnan(myTrajectory.centerline), 2), :);

plot(myTrajectory.centerline(:, 1), myTrajectory.centerline(:, 2), 'r-', 'LineWidth', 2);

plot(centerline(:, 1), centerline(:, 2), 'r-', 'LineWidth', 2);

num_points = size(myTrajectory.centerline, 1);

%text(myTrajectory.centerline(1, 1), myTrajectory.centerline(1, 2), ['Number of Points: ', num2str(num_points)], 'Color', 'g');

grid on;

% Figure configurations

set(gca, 'fontsize', FontSize_axis, 'FontName', 'Times', 'fontweight', 'bold')

title('\bf{Tracking Trajectory}', ...

'FontUnits', 'points', ...

'FontWeight', font_weight, ...

'interpreter', 'latex', ...

'FontSize', FontSize, ...

'FontName', 'Times')

xlabel('\bf{x(m)}', ...

'FontUnits', 'points', ...

'FontWeight', font_weight, ...

'interpreter', 'latex', ...

'FontSize', FontSize, ...

'FontName', 'Times')

ylabel('\bf{y(m)}', ...

'FontUnits', 'points', ...

'FontWeight', font_weight, ...

'interpreter', 'latex', ...

'FontSize', FontSize, ...

'FontName', 'Times')

% Make the gif file

frame = getframe(h);

img = frame2im(frame);

[AA, map] = rgb2ind(img, 256);

if i == 1

imwrite(AA, map, filename, 'gif', 'LoopCount', Inf, 'DelayTime', delay_time);

else

imwrite(AA, map, filename, 'gif', 'WriteMode', 'append', 'DelayTime', delay_time);

end

% if ~any(isnan([min_x, max_x, min_y, max_y])) & ~any(isinf([min_x, max_x, min_y, max_y]))

% % Entferne NaN-Werte aus der Centerline

% centerline_clean = centerline(~any(isnan(centerline), 2), :);

% Finden Sie die Minimum- und Maximumwerte von centerline

if ~any(isnan([min_x, max_x, min_y, max_y])) & ~any(isinf([min_x, max_x, min_y, max_y]))

% Überprüfen Sie, ob die Werte für min und max konsistent sind

centerline_clean = centerline(~any(isnan(centerline), 2), :);

min_x = min(centerline(:, 1), 'omitnan'); % Ignoriere NaN-Werte

min_y = min(centerline(:, 2), 'omitnan'); % Ignoriere NaN-Werte

max_x = max(centerline(:, 1), 'omitnan'); % Ignoriere NaN-Werte

max_y = max(centerline(:, 2), 'omitnan'); % Ignoriere NaN-Werte

tolerance = 1e-9; % Toleranzwert für den Mindestabstand

if min_x + tolerance < max_x & min_y + tolerance < max_y

% Setze xlim und ylim

xlim([min_x - win, max_x + win]);

ylim([min_y - win, max_y + win]);

else

disp('Warnung: Die Centerline enthält möglicherweise ungültige Werte für Achsenlimits.');

disp(sprintf('min_x: %f, max_x: %f', min_x, max_x));

disp(sprintf('min_y: %f, max_y: %f', min_y, max_y));

end

else

disp('Warnung: Die Centerline enthält ungültige Werte für Achsenlimits.');

disp(['min_x: ', num2str(min_x), ', max_x: ', num2str(max_x)]);

disp(['min_y: ', num2str(min_y), ', max_y: ', num2str(max_y)]);

end

[Ivan Bullier Ananfack Tsakem]

classdef Car<handle

properties

% States of the Model

states; % [x,y,phi,velocity]

augmented_states;

control_inputs; %[steering_angle, acceleration]

% Geometry of the car

track=2057/10000; %meter

wheel_base=2665/10000; %meter

wheel_diameter=66.294/1000;

wheel_width=25.908/1000;

%Time step for the simulation

ts=0.1;

% Model limits

steering_angle_limit = deg2rad(33.75);

max_velocity = 98.3488; % m/s

% PID Errors

old_cte;

cte_integral;

end

methods

%% Constructor

function obj=Car(states,control_inputs)

% Check not to exceed the maximum velocity.

if (states(4) > obj.max_velocity)

states(4)= obj.max_velocity;

end

obj.states=states;

% Normalize the steering angle to be between -pi and +pi.

control_inputs(1) = atan2(sin(control_inputs(1)),cos(control_inputs(1)));

% Check that the steering angle is not exceeding the limit.

if (control_inputs(1) > obj.steering_angle_limit)

control_inputs(1) = obj.steering_angle_limit;

elseif (control_inputs(1) < -obj.steering_angle_limit)

control_inputs(1) = -obj.steering_angle_limit;

end

obj.control_inputs=control_inputs;

obj.old_cte = 0;

obj.cte_integral = 0;

end

%% Plot My car

function show(obj)

%Plot the Body of the Car

pc1=[-obj.wheel_base/2;obj.track/2;1];

pc2=[obj.wheel_base/2;obj.track/2;1];

pc3=[obj.wheel_base/2;-obj.track/2;1];

pc4=[-obj.wheel_base/2;-obj.track/2;1];

pc=[pc1 pc2 pc3 pc4];

pw1=[-obj.wheel_diameter/2;obj.wheel_width/2;1];

pw2=[obj.wheel_diameter/2;obj.wheel_width/2;1];

pw3=[obj.wheel_diameter/2;-obj.wheel_width/2;1];

pw4=[-obj.wheel_diameter/2;-obj.wheel_width/2;1];

pwheel=[pw1 pw2 pw3 pw4];

[x,y,~,~]=state_unpack(obj);

%Plot the center of the car

plot(x,y,'o','Markersize',10,'Markerface','b')

hold on

R_car_world=carf_to_wf(obj);

pcw=R_car_world*pc;

plot([pcw(1,1) pcw(1,2)],[pcw(2,1) pcw(2,2)],'b','Linewidth',2)

plot([pcw(1,2) pcw(1,3)],[pcw(2,2) pcw(2,3)],'b','Linewidth',2)

plot([pcw(1,3) pcw(1,4)],[pcw(2,3) pcw(2,4)],'b','Linewidth',2)

plot([pcw(1,4) pcw(1,1)],[pcw(2,4) pcw(2,1)],'b','Linewidth',2)

%plot the Back wheels

%Plot the Front wheels

[si,~]=control_inputs_unpack(obj);

for i=2:3

R_wheel_to_car=wheel_frame_to_car(pc(:,i),si);

pwheel=R_car_world*R_wheel_to_car*pwheel;

plot([pwheel(1,1) pwheel(1,2)],[pwheel(2,1) pwheel(2,2)],'Color','red','Linewidth',2)

plot([pwheel(1,2) pwheel(1,3)],[pwheel(2,2) pwheel(2,3)],'Color','red','Linewidth',2)

plot([pwheel(1,3) pwheel(1,4)],[pwheel(2,3) pwheel(2,4)],'Color','red','Linewidth',2)

plot([pwheel(1,4) pwheel(1,1)],[pwheel(2,4) pwheel(2,1)],'Color','red','Linewidth',2)

pwheel=[pw1 pw2 pw3 pw4];

end

for i=1:3:4

R_wheel_to_car=wheel_frame_to_car(pc(:,i),0);

pwheel=R_car_world*R_wheel_to_car*pwheel;

plot([pwheel(1,1) pwheel(1,2)],[pwheel(2,1) pwheel(2,2)],'Color','red','Linewidth',2)

plot([pwheel(1,2) pwheel(1,3)],[pwheel(2,2) pwheel(2,3)],'Color','red','Linewidth',2)

plot([pwheel(1,3) pwheel(1,4)],[pwheel(2,3) pwheel(2,4)],'Color','red','Linewidth',2)

plot([pwheel(1,4) pwheel(1,1)],[pwheel(2,4) pwheel(2,1)],'Color','red','Linewidth',2)

pwheel=[pw1 pw2 pw3 pw4];

end

set(gca,'units','centimeters')

hold off

end

%% Get the transfomration matrix from car coordinate system to world coordinate system.

function R=carf_to_wf(obj)

[x,y,phi,~]=state_unpack(obj);

R=[cos(phi) -sin(phi) x;sin(phi) cos(phi) y;0 0 1];

end

%% Unpack the states

function [x,y,phi,v]=state_unpack(obj)

x=obj.states(1);

y=obj.states(2);

phi=obj.states(3);

v=obj.states(4);

end

%% Unpack control_inputs

function [si,acc]=control_inputs_unpack(obj)

si=obj.control_inputs(1);

acc=obj.control_inputs(2);

end

%% The dynamics of the car

% This function predicts the next state of the car based on its

% current state and the control inputs

function obj=update_state(obj)

[x,y,phi,v] = state_unpack(obj);

[si,acc] = control_inputs_unpack(obj);

x_next = x + v*cos(phi)*obj.ts;

y_next = y + v*sin(phi)*obj.ts;

phi_next = phi + v/(obj.wheel_base)*si*obj.ts;

v_next = v + acc*obj.ts;

% Check not the exceed the speed limit.

if(v_next > obj.max_velocity)

v_next = obj.max_velocity;

end

obj.states=[x_next y_next phi_next v_next];

end

%% A simple setter function. It is not necessary, but it will make the main code more readable.

function obj = update_input(obj, control_inputs)

% Normalize the steering angle to be between -pi and +pi.

control_inputs(1) = atan2(sin(control_inputs(1)), cos(control_inputs(1)));

% Check that the steering angle is not exceeding the limit.

if (control_inputs(1) > obj.steering_angle_limit)

control_inputs(1) = obj.steering_angle_limit;

elseif (control_inputs(1) < -obj.steering_angle_limit)

control_inputs(1) = -obj.steering_angle_limit;

end

obj.control_inputs = control_inputs;

end

function PID_Controller(obj, cte)

kp = 0.02;

kd = 0.25;

ki = 0.00005;

dcte = cte - obj.old_cte;

obj.cte_integral = obj.cte_integral + cte;

obj.old_cte = cte;

steering = kp * cte + kd * dcte + ki * obj.cte_integral;

[~, a] = obj.control_inputs_unpack;

control_signal = [steering, a];

obj.update_input(control_signal);

end

end

end

%% Helper Functions

function R=wheel_frame_to_car(pc,si)

R=[cos(si) -sin(si) pc(1);sin(si) cos(si) pc(2);0 0 1];

end

[Walter Roberson]

 car_states = Car.state_unpack();

but state_unpack is defined as returning up to four parts, starting with x and then y.

You then proceed to index car_states, which is equivalent to indexing the returned x. Do we have solid reason to expect that the x will always have at least two components?

[Ivan Bullier Ananfack Tsakem]

bitte ich hätte gerne diesen PID-Regler für alle Stecke implimentieren und optimieren. können Sie mir bitte helfen?

classdef Car<handle

properties

% States of the Model

states; % [x,y,phi,velocity]

augmented_states;

control_inputs; %[steering_angle, acceleration]

% Geometry of the car

track=2057/1000; %meter

wheel_base=2665/1000; %meter

wheel_diameter=66.294/100;

wheel_width=25.908/100;

%Time step for the simulation

ts=0.05;

% Model limits

%steering_angle_limit = deg2rad(33.75);

steering_angle_limit = deg2rad(60.75);

max_velocity = 40.3488; % m/s

%max_velocity = 98.3488; % m/s

% PID Errors

old_cte;

cte_intergral;

last_smoothed_inputs = [0 0]

max_acceleration = 6.0;

end

methods

%% Constructor

function obj=Car(states,control_inputs)

% Check not to exceed the maximum velocity.

if (states(4) > obj.max_velocity)

states(4)= obj.max_velocity;

end

obj.states=states;

% Normalize the steering angle to be between -pi and +pi.

control_inputs(1) = atan2(sin(control_inputs(1)),cos(control_inputs(1)));

% Check that the steering angle is not exceeding the limit.

if (control_inputs(1) > obj.steering_angle_limit)

control_inputs(1) = obj.steering_angle_limit;

elseif (control_inputs(1) < -obj.steering_angle_limit)

control_inputs(1) = -obj.steering_angle_limit;

end

obj.control_inputs=control_inputs;

obj.old_cte = 0;

obj.cte_intergral = 0;

end

%% Plot My car

function show(obj)

%Plot the Body of the Car

pc1=[-obj.wheel_base/2;obj.track/2;1];

pc2=[obj.wheel_base/2;obj.track/2;1];

pc3=[obj.wheel_base/2;-obj.track/2;1];

pc4=[-obj.wheel_base/2;-obj.track/2;1];

pc=[pc1 pc2 pc3 pc4];

pw1=[-obj.wheel_diameter/2;obj.wheel_width/2;1];

pw2=[obj.wheel_diameter/2;obj.wheel_width/2;1];

pw3=[obj.wheel_diameter/2;-obj.wheel_width/2;1];

pw4=[-obj.wheel_diameter/2;-obj.wheel_width/2;1];

pwheel=[pw1 pw2 pw3 pw4];

[x,y,~,~]=state_unpack(obj);

%Plot the center of the car

plot(x,y,'o','Markersize',10,'Markerface','b')

hold on

R_car_world=carf_to_wf(obj);

pcw=R_car_world*pc;

plot([pcw(1,1) pcw(1,2)],[pcw(2,1) pcw(2,2)],'b','Linewidth',2)

plot([pcw(1,2) pcw(1,3)],[pcw(2,2) pcw(2,3)],'b','Linewidth',2)

plot([pcw(1,3) pcw(1,4)],[pcw(2,3) pcw(2,4)],'b','Linewidth',2)

plot([pcw(1,4) pcw(1,1)],[pcw(2,4) pcw(2,1)],'b','Linewidth',2)

%plot the Back wheels

%Plot the Front wheels

[si,~]=control_inputs_unpack(obj);

for i=2:3

R_wheel_to_car=wheel_frame_to_car(pc(:,i),si);

pwheel=R_car_world*R_wheel_to_car*pwheel;

plot([pwheel(1,1) pwheel(1,2)],[pwheel(2,1) pwheel(2,2)],'Color','red','Linewidth',2)

plot([pwheel(1,2) pwheel(1,3)],[pwheel(2,2) pwheel(2,3)],'Color','red','Linewidth',2)

plot([pwheel(1,3) pwheel(1,4)],[pwheel(2,3) pwheel(2,4)],'Color','red','Linewidth',2)

plot([pwheel(1,4) pwheel(1,1)],[pwheel(2,4) pwheel(2,1)],'Color','red','Linewidth',2)

pwheel=[pw1 pw2 pw3 pw4];

end

for i=1:3:4

R_wheel_to_car=wheel_frame_to_car(pc(:,i),0);

pwheel=R_car_world*R_wheel_to_car*pwheel;

plot([pwheel(1,1) pwheel(1,2)],[pwheel(2,1) pwheel(2,2)],'Color','red','Linewidth',2)

plot([pwheel(1,2) pwheel(1,3)],[pwheel(2,2) pwheel(2,3)],'Color','red','Linewidth',2)

plot([pwheel(1,3) pwheel(1,4)],[pwheel(2,3) pwheel(2,4)],'Color','red','Linewidth',2)

plot([pwheel(1,4) pwheel(1,1)],[pwheel(2,4) pwheel(2,1)],'Color','red','Linewidth',2)

pwheel=[pw1 pw2 pw3 pw4];

end

set(gca,'units','centimeters')

hold off

end

%% Get the transfomration matrix from car coordinate system to world coordinate system.

function R=carf_to_wf(obj)

[x,y,phi,~]=state_unpack(obj);

R=[cos(phi) -sin(phi) x;sin(phi) cos(phi) y;0 0 1];

end

%% Unpack the states

function [x,y,phi,v]=state_unpack(obj)

x=obj.states(1);

y=obj.states(2);

phi=obj.states(3);

v=obj.states(4);

end

%% Unpack control_inputs

function [si,acc]=control_inputs_unpack(obj)

si=obj.control_inputs(1);

acc=obj.control_inputs(2);

end

%% The dynamics of the car

% This function predicts the next state of the car based on its

% current state and the control inputs

function obj=update_state(obj)

[x,y,phi,v] = state_unpack(obj);

[si,acc] = control_inputs_unpack(obj);

x_next = x + v*cos(phi)*obj.ts;

y_next = y + v*sin(phi)*obj.ts;

phi_next = phi + v/(obj.wheel_base)*si*obj.ts;

v_next = v + acc*obj.ts;

% Check not the exceed the speed limit.

if(v_next > obj.max_velocity)

v_next = obj.max_velocity;

end

obj.states=[x_next y_next phi_next v_next];

end

%% A simple setter function. It is not necessary, but it will make the main code more readable.

function obj=update_input(obj,control_inputs)

% Normalize the steering angle to be between -pi and +pi.

control_inputs(1) = atan2(sin(control_inputs(1)),cos(control_inputs(1)));

% Check that the steering angle is not exceeding the limit.

if (control_inputs(1) > obj.steering_angle_limit)

control_inputs(1) = obj.steering_angle_limit;

elseif (control_inputs(1) < -obj.steering_angle_limit)

control_inputs(1) = -obj.steering_angle_limit;

end

obj.control_inputs=control_inputs;

end

function smooth_control_inputs(obj, smoothing_factor)

% Smooth control inputs using a simple low-pass filter

current_inputs = obj.control_inputs;

smoothed_inputs = smoothing_factor * current_inputs + (1 - smoothing_factor) * obj.last_smoothed_inputs;

obj.last_smoothed_inputs = smoothed_inputs;

obj.control_inputs = smoothed_inputs;

end

% function LateralLongitudinal_Controller(obj, cte, desired_yaw_rate)

% kp_lat = 0.03; % Parameter für laterale Regelung

% kd_lat = 0.3;

% ki_lat = 0.0001;

%

% kp_lon = 0.1; % Parameter für longitudinale Regelung

% kd_lon = 0.01;

% ki_lon = 0.0001;

%

% dcte = cte - obj.old_cte;

% obj.cte_intergral = obj.cte_intergral + cte;

% obj.old_cte = cte;

%

% % Lateral Control (quer)

% steering = kp_lat * cte + kd_lat * dcte + ki_lat * obj.cte_intergral;

%

% % Longitudinal Control (längs)

% desired_velocity = obj.max_velocity; % Setzen Sie die gewünschte Geschwindigkeit

% velocity_error = desired_velocity - obj.states(4);

%

% % Fügen Sie die max_acceleration-Eigenschaft hinzu

% max_acceleration = 2.0; % Setzen Sie hier den gewünschten maximalen Beschleunigungswert ein

%

% acceleration = kp_lon * velocity_error + kd_lon * obj.states(4) + ki_lon * sum(obj.states(1:4));

%

% % Begrenzen Sie die Beschleunigung

% acceleration = max(-max_acceleration, min(max_acceleration, acceleration));

%

% control_signal = [steering, acceleration];

%

% obj.update_input(control_signal);

% end

function PID_Controller(obj,cte)

kp = 0.08;

kd = 0.4;

ki = 0.000001;

dcte = cte - obj.old_cte;

obj.cte_intergral = obj.cte_intergral + cte;

obj.old_cte = cte;

steering = kp * cte + kd * dcte + ki * obj.cte_intergral;

[~,a] = obj.control_inputs_unpack;

control_signal = [steering,a];

obj.update_input(control_signal);

end

end

end

%% Helper Functions

function R=wheel_frame_to_car(pc,si)

R=[cos(si) -sin(si) pc(1);sin(si) cos(si) pc(2);0 0 1];

end

hier ist die trajectory classe:

classdef Trajectory < handle

properties

centerline;

nump = 8;

Nearst_Points;

Nearst_Points_C;

cte;

show_nearest = true;

show_fit = true;

Coeff_fit;

end

methods

function obj = Trajectory(file_path)

% Read Pylon data from file

pylons_data = readtable(file_path, 'Delimiter', ',', 'Format', '%s%f%f');

blue_pylons = pylons_data(strcmp(pylons_data.Var1, 'blue'), :);

yellow_pylons = pylons_data(strcmp(pylons_data.Var2, 'yellow'), :);

% Use pylons_data to generate centerline (replace this with your logic)

obj.centerline = generate_trajectory(pylons_data);

obj.Nearst_Points = zeros(obj.nump, 2);

obj.Nearst_Points_C = zeros(obj.nump, 2);

obj.cte = 0;

obj.Coeff_fit = [0 0 0];

end

function show(obj, Car)

% No need to read pylons_data from file since you're using centerline

% Überprüfen Sie, ob die Centerline leer oder gültig ist

if ~isempty(obj.centerline) && size(obj.centerline, 1) > 1

% Plot Centerline

plot(obj.centerline(:, 1), obj.centerline(:, 2), 'LineWidth', 1.5);

hold on;

else

% Die Centerline ist leer oder ungültig, geben Sie eine Warnung aus

warning('Die Centerline ist leer oder nicht korrekt zugewiesen.');

end

% Hier den Zustand des Autos abrufen und die Variable carPositionOnCenterline aktualisieren

Car_states = Car.state_unpack();

% Deine weiteren Berechnungen und Plots hier ...

hold off;

end

function W_to_Car_Coordinate_system(obj, Car)

[x, y, phi, ~] = Car.state_unpack;

tranlate = obj.Nearst_Points - [x, y];

obj.Nearst_Points_C = ([cos(phi) sin(phi); -sin(phi) cos(phi)] * tranlate')';

end

function find_nearest_points(obj, Car)

[x, y, ~, ~] = Car.state_unpack;

dist_list = sqrt((obj.centerline(:, 1) - x).^2 + (obj.centerline(:, 2) - y).^2);

[~, indices] = sort(dist_list);

obj.Nearst_Points = obj.centerline(indices(1:obj.nump), :);

end

function poly_fit(obj, Car)

obj.W_to_Car_Coordinate_system(Car);

obj.Coeff_fit = polyfit(obj.Nearst_Points_C(:, 1), obj.Nearst_Points_C(:, 2), 2);

end

function compute_error(obj)

obj.cte = obj.Coeff_fit(end);

end

end

end

% Transform points to car coordinate system

function CarPoints = CarToCarCoordinateSystem(WorldPoints, car)

[x, y, phi, ~] = car.state_unpack;

rotate = ([cos(phi), -sin(phi); sin(phi), cos(phi)] * WorldPoints')';

CarPoints = rotate + [x, y];

end

so dunktioniert der code: