Index exceeds number of array elements (1)?
Mostrar comentarios más antiguos
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?
1 comentario
Mark Linne
el 30 de Jul. de 2021
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)?
Respuesta aceptada
Más respuestas (9)
sura naji
el 18 de Feb. de 2020
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)]?
1 comentario
BAMIDELE JEGEDE
el 31 de Mayo de 2020
Editada: Walter Roberson
el 4 de Jun. de 2021
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
3 comentarios
Aquatris
el 31 de Mayo de 2020
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?
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
el 12 de Feb. de 2021
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
Sara Babaie
el 12 de Feb. de 2021
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)
3 comentarios
Walter Roberson
el 12 de Feb. de 2021
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
el 12 de Feb. de 2021
This error is related to the dimensions and the number of elements is more.
Walter Roberson
el 12 de Feb. de 2021
In this case, Sara's error is in forgetting a multiplication symbol. stfb*(t84^4-tsurr^4)
Simon Hudon-Labonté
el 9 de Mzo. de 2021
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
4 comentarios
Farshid R
el 9 de Mzo. de 2021
This error is related to the dimensions and the number of elements is more.
Walter Roberson
el 9 de Mzo. de 2021
You would get that error if the first parameter was an array instead of a function handle
Binyam Sime
el 3 de Jun. de 2021
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
el 3 de Jun. de 2021
sum =sum +freq(i);
You do not define freq() in the code.
Wu Changjin Wu
el 4 de Jun. de 2021
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)
1 comentario
Walter Roberson
el 4 de Jun. de 2021
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)
ELIZABETH ESPARZA
el 5 de Jun. de 2021
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
5 comentarios
Walter Roberson
el 5 de Jun. de 2021
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
el 5 de Jun. de 2021
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
el 5 de Jun. de 2021

ELIZABETH ESPARZA
el 5 de Jun. de 2021
oh well, but the mistake is in line 58
Aquatris
el 6 de Jun. de 2021
That index problem happens because of variable w. You never assigned values to w(i+1) in your code.
emre bahçeci
el 15 de Jun. de 2021
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
2 comentarios
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]
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,1)=subs(Q(x),{theta_dot_dot theta_dot theta},{tetai(i) tetah(i) teta(i)});
end
Q_new
emre bahçeci
el 16 de Jun. de 2021
Thanks a lot
Anastasiya Moiseeva
el 23 de Mzo. de 2023
Editada: Anastasiya Moiseeva
el 23 de Mzo. de 2023
0 votos
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
%----------------------------------------
1 comentario
Walter Roberson
el 23 de Mzo. de 2023
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.
labtst
el 3 de En. de 2024
0 votos
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)
>>
5 comentarios
labtst
el 3 de En. de 2024
Movida: Walter Roberson
el 3 de En. de 2024
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
%
labtst
el 3 de En. de 2024
Movida: Walter Roberson
el 3 de En. de 2024
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
labtst
el 3 de En. de 2024
Movida: Walter Roberson
el 3 de En. de 2024
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
el 3 de En. de 2024
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?
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:
Categorías
Más información sobre Axis Labels en Centro de ayuda y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
