Failure in initial nonlinear constraint function evaluation. FMINCON cannot continue.
Show older comments
The following is the script which I run the fmincon solver in. I am not sure whether the error is resulting because the non-linear constraint is incorrectly defined or whether my initial value x0 does not satisfy the constraints.
X = [0,0,0];
u = [0,0,0];
v_fb = [-4.2426,0,-0.7854];
n = 9;
deltaT = 0.1;
n_t = 1;
X_h_o = [2,2,0,0];
pos = safetyModule(X, u, v_fb, n, deltaT, n_t, X_h_o)
function pos = safetyModule(X, u, v_fb, n, deltaT, n_t, X_h_o)
%extend system state to match number of discretized time points in
%horizon including t = 0
X_h = [X, u, zeros(1, 6 * (n-1))];
%initial guess (may not be feasible)
x0 = X_h;
%define constants
D_obs = 0.5; %margin of safety
alpha = 1;
%adjust variables
n_t = n_t + 1;
%decompose v_fb
xdot_fb = v_fb(1);
ydot_fb = v_fb(2);
w_fb = v_fb(3);
%objective function
fun = @(X_h) 0;
for i = 1:n
fun1 = abs(X_h(i+3*n) - xdot_fb) + abs(X_h(i+4*n) - ydot_fb) + abs(X_h(i+5*n) - w_fb);
fun = @(X_h) fun(X_h) + fun1;
end
%equality constraints
Aeq = zeros(3*n, 6 * n);
Aeq(1,1) = 1;
Aeq(2,1+n) = 1;
Aeq(3,1+2*n) = 1;
for i = 1:n-1
Aeq(3+i, i) = -1; Aeq(3+i,i+1) = 1; Aeq(3+i,i+n*3) = -deltaT;% x dynamics
Aeq(2 + n + i, n + i) = -1; Aeq(2 + n + i, n + i + 1) = 1; Aeq(2 + n + i, 4 * n + i) = -deltaT;% y dynamics
Aeq(1 + 2*n + i, 2* n + i) = -1; Aeq(1 + 2*n + i, 2 * n + i + 1) = 1; Aeq(1 + 2*n + i, 5 * n + i) = -deltaT;% z dynamics
end
beq = zeros(3*n, 1);
beq(1)=X_h(1);
beq(2)=X_h(2);
beq(3)=X_h(3);
%inequality constraints
nonlcon = @(X_h) CBFcon(X_h, n, deltaT, X_h_o, alpha, D_obs);
A = [];
b = [];
%bounds
ub = [inf(3*n, 1), 2.5 * ones(3*n, 1)];
lb = [-inf(3*n, 1), -2.5 * ones(3*n, 1)];
%run the solver
dec = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon);
%extract xdot, ydot, and w for first n_p timesteps
pos = zeros(1, (n_t - 1)* 3);
j = 1;
for i = 2:n_t
pos(j) = dec(i);
pos(j+1) = dec(i+n);
pos(j+2) = dec(i+2*n);
j = j + 3;
end
end
I define my nonlinear constraint in a seperate file:
function [c, ceq] = CBFcon(X, n, deltaT, X_o, alpha, D_obs)
x_curr = X_o(1);
y_curr = X_o(3);
x_vel = X_o(2);
y_vel = X_o(4);
for i = 1:n
c(i) = -(((X(i) - x_curr)/sqrt((X(i) - x_curr)^2 + (X(i+n) - y_curr)^2 )) * X(i + 3*n) + ((X(i+n) - y_curr)/sqrt((X(i) - x_curr)^2 + (X(i+n) - y_curr)^2 )) * X(i + 4*n) + alpha(sqrt((X(i) - x_curr)^2 + (X(i+n) - y_curr)^2 ) - D_obs));
x_curr = x_curr + x_vel * deltaT;
y_curr = y_curr + y_vel * deltaT;
end
ceq = [];
end
4 Comments
Torsten
on 23 Mar 2024
I am not sure whether the error is resulting because the non-linear constraint is incorrectly defined or whether my initial value x0 does not satisfy the constraints.
We don't know, either, because you didn't supply the script part in which you call "safetyModule" with appropriate input values.
I am not sure whether the error is resulting because the non-linear constraint is incorrectly defined or whether my initial value x0 does not satisfy the constraints.
The former. But why not just check it yourself by calling your constraint function,
[c,ceq]=nonlcon(x0)
Dhruv
on 23 Mar 2024
Yes, but providing the values as a screenshot will not work because -
- It cannot be copy-pasted.
- It doesn't gives us the values accurately, because the screenshot only displays them to four decimal places
Please attach a .mat file with the values.
Answers (1)
Torsten
on 23 Mar 2024
0 votes
Most probably you mean
c(i) = -(((X(i) - x_curr)/sqrt((X(i) - x_curr)^2 + (X(i+n) - y_curr)^2 )) * X(i + 3*n) + ((X(i+n) - y_curr)/sqrt((X(i) - x_curr)^2 + (X(i+n) - y_curr)^2 )) * X(i + 4*n) + alpha*(sqrt((X(i) - x_curr)^2 + (X(i+n) - y_curr)^2 ) - D_obs));
instead of
c(i) = -(((X(i) - x_curr)/sqrt((X(i) - x_curr)^2 + (X(i+n) - y_curr)^2 )) * X(i + 3*n) + ((X(i+n) - y_curr)/sqrt((X(i) - x_curr)^2 + (X(i+n) - y_curr)^2 )) * X(i + 4*n) + alpha(sqrt((X(i) - x_curr)^2 + (X(i+n) - y_curr)^2 ) - D_obs));
25 Comments
Dhruv
on 23 Mar 2024
Gradient decent is the simplest possible algorithm for optimization. It's not implemented in "fmincon".
Dhruv
on 23 Mar 2024
Dhruv
on 23 Mar 2024
Torsten
on 23 Mar 2024
for i = 1:n
fun1 = @(X_h) abs(X_h(i+3*n) - xdot_fb) + abs(X_h(i+4*n) - ydot_fb) + abs(X_h(i+5*n) - w_fb);
fun = @(X_h) fun(X_h) + fun1(X_h);
end
instead of
for i = 1:n
fun1 = abs(X_h(i+3*n) - xdot_fb) + abs(X_h(i+4*n) - ydot_fb) + abs(X_h(i+5*n) - w_fb);
fun = @(X_h) fun(X_h) + fun1;
end
And if you mean variables that can only take discrete values by "decision variables", then you cannot solve this problem using "fmincon". "fmincon" can only solve optimization problems where all variables are continuous.
Torsten
on 24 Mar 2024
The solution strategy of "fmincon" is based on the differentiability of objective and constraints with respect to the solution variables. This is not given if some or all of the solution variables are discrete.
You will have to use "ga" instead with the discrete variables defined in the "intcon" array. Nonlinear constraints are only allowed as inequalities, but this seems to be the case in your problem.
Dhruv
on 26 Mar 2024
Dhruv
on 26 Mar 2024
Dhruv
on 26 Mar 2024
Dhruv
on 26 Mar 2024
In the code below, equality and inequality constraints imposed by Aeq, beq and A,b are respected.
Please give a counterexample.
X = [0,0]; %starting state
u = [0,0];%input at t = 0
%starting state of obstacle:
X_o = zeros(1,4);
X_o(1) = 0;
X_o(2) = 4;
% X_o(3) = 0;
% X_o(4) = -0.4;
position_accuracy = 0.05;
X_des = [3,3]; %desired position
t_max = 3; %maximum simulation time
n = 2; %horizon
deltaT = 0.05;%timestem length
t = 1; %starting time
goal(1,1) = 10;
goal(1,2) = 10;
while norm(goal(1,2) - X(1,2)) > position_accuracy && t <= t_max
%describe behavior of obstacle
% if t > 10 && t <= 12
% X_o(3) = -2;
% X_o(4) = 0;
% end
%
% if t > 12
% X_o(3) = 1;
% X_o(4) = -1;
% end
v_fb = trajectoryController_h(X, X_des);
%fprintf('TIME STEP: %d\n', t);
[X, u] = safetyModule_h(X, u, v_fb, n, deltaT, X_o);
%can use input to determine theta
theta = atan(u(2)/u(1));
%print metrics at each timestep:
% fprintf('Input: ');
% for i = 1:length(u)
% fprintf('%d', u(i));
% end
% fprintf('\n');
%
% fprintf('Feedback velocities: ');
% for i = 1:length(v_fb)
% fprintf('%d ', v_fb(i));
% end
% fprintf('\n');
%
% fprintf('Robot State: ');
% for i = 1:length(X)
% fprintf('%d ', X(i));
% end
% fprintf('\n \n');
%archive and plot
X_graph(t) = X(1);
Y_graph(t) = X(2);
cla;
daspect([1 1 1]);
hold on;
plot(X_o(1), X_o(2), 'or', 'MarkerSize', 5);
plot(goal(1), goal(2), 'og', 'MarkerSize', 5);
plot(X_graph(1:t), Y_graph(1:t), '-b'); %traveled path
plot([X(1) X(1)+0.2*cos(theta)], [X(2) X(2)+0.2*sin(theta)], '-r'); %heading of robot
xlim([0, 4]);
ylim([0, 4]);
drawnow;
pause(deltaT);
%iterate
t = t + 1;
X_o(1) = X_o(1) + X_o(3)* deltaT;
X_o(2) = X_o(2) + X_o(4)* deltaT;
end
%{
trackingController calculates the desired velocity based on current,
position, desired position, and previous timestep input
%}
function v_fb = trajectoryController_h(X, X_des)
%define proportional constants for P controller
K_x = 1;
K_y = 1;
x = X(1);
y = X(2);
x_des = X_des(1);
y_des = X_des(2);
%extract errors
ex = x_des - x;
ey = y_des - y;
%Proportional controller
vfb_x = K_x * ex;
vfb_y = K_y * ey;
v_fb = [vfb_x vfb_y];
%final stretch constant speed definition:
if (norm(X_des - x,y) < 0.75)
mag = norm(v_fb);
v_fb = (v_fb/mag) * 1;
end
end
function [pos, u] = safetyModule_h(X, u, v_fb, n, deltaT, X_h_o)
%extend system state to match number of discretized time points in
%horizon including t = 0
X_h = [X, u, ones(1, 4 * (n-1))];
% fprintf('X_h: %d\n', X_h(1:4));
%initial guess (may not be feasible)
x0 = X_h;
%define constants
D_obs = 0.1; %margin of safety
alpha = 1;
%decompose v_fb
xdot_fb = v_fb(1);
ydot_fb = v_fb(2);
%objective function
fun = @(X_h) 0;
for i = 1:n
fun1 = @(X_h) abs(X_h(i+2*n) - xdot_fb) + abs(X_h(i+3*n) - ydot_fb);
fun = @(X_h) fun(X_h) + fun1(X_h);
end
%equality constraints
Aeq = zeros(2 * n, 4 * n);
Aeq(1,1) = 1;
Aeq(2,1+n) = 1;
for i = 1:n-1
Aeq(2+i, i) = -1; Aeq(2+i,i+1) = 1; Aeq(2+i,i+n*2) = -deltaT;% x dynamics
Aeq(1 + n + i, n + i) = -1; Aeq(1 + n + i, n + i + 1) = 1; Aeq(1 + n + i, 3 * n + i) = -deltaT;% y dynamics
end
beq = zeros(2*n, 1);
beq(1)=X_h(1);
beq(2)=X_h(2);
% fprintf('Aeq:\n');
% for i = 1:size(Aeq, 1)
% for j = 1:size(Aeq, 2)
% fprintf('%d\t', Aeq(i, j)); % Change '%d' to '%f' for floating-point numbers
% end
% fprintf('\n'); % New line for each row
% end
%inequality constraints
nonlcon = @(X_h) CBFcon_h(X_h, n, deltaT, X_h_o, alpha, D_obs);
A = zeros(4*(n-1),4*n);
for i = 1:n-1
%less than or equal to acceleration constraint
A(i, 2*n + i) = -1; A(i, 2*n + i + 1) = 1;
A((n-1) + i, 3*n + i) = -1; A((n-1) + i, 3*n + i + 1) = 1;
%greater than or equal to acceleration
A(2*(n-1) + i, 2*n + i) = 1; A(2*(n-1) + i, 2*n + i + 1) = -1;
A(3*(n-1) + i, 3*n + i) = 1; A(3*(n-1) + i, 3*n + i + 1) = -1;
end
b = deltaT * 0.5 * ones(4*(n-1), 1);
% fprintf('A:\n');
% for i = 1:size(A, 1)
% for j = 1:size(A, 2)
% fprintf('%d\t', A(i, j)); % Change '%d' to '%f' for floating-point numbers
% end
% fprintf('\n'); % New line for each row
% end
%bounds
ub = [4 * ones(2*n, 1), 2.5 * ones(2*n, 1)];
lb = [zeros(2*n, 1), -2.5 * ones(2*n, 1)];
%run the solver
options = optimoptions('fmincon', 'Display', 'off');
dec = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon, options);
A*dec.'-b
Aeq*dec.'-beq
%extract x, y for the upcoming timestep
pos = zeros(1, 2);
pos(1) = dec(2);
pos(2) = dec(2+n);
%extract xdot, ydot for upcoming timestep
u = zeros(1, 2);
u(1) = dec(2*n + 2);
u(2) = dec(3*n + 2);
% fprintf('xdot_fb: %d\n', xdot_fb);
% fprintf('ydot_fb: %d\n', ydot_fb);
% fprintf('u(1): %d\n', u(1));
% fprintf('u(2): %d\n', u(2));
end
function [c, ceq] = CBFcon_h(X_h, n, deltaT, X_o, alpha, D_obs)
x_curr = X_o(1);
y_curr = X_o(2);
x_vel = X_o(3);
y_vel = X_o(4);
c = zeros(1, n);
% Loop through iterations
for i = 1:n
% Calculate distance between X_h(i) and current position
d = sqrt((X_h(i) - x_curr)^2 + (X_h(i + n) - y_curr)^2);
% Calculate A and B
A = (X_h(i) - x_curr) / d;
B = (X_h(i + n) - y_curr) / d;
% Calculate constraint value for c(i)
c(i) = - (A * X_h(i + 2*n) + B * X_h(i + 3*n) + alpha * (d - D_obs));
% Update current positions using velocity and time step
x_curr = x_curr + x_vel * deltaT;
y_curr = y_curr + y_vel * deltaT;
end
% Return empty equality constraint array
ceq = [];
end
Dhruv
on 27 Mar 2024
Changing the solution method seems to work here:
X = [0,0]; %starting state
u = [0,0];%input at t = 0
%starting state of obstacle:
X_o = zeros(1,4);
X_o(1) = 2;
X_o(2) = 2;
X_o(3) = 0;
X_o(4) = 0;
position_accuracy = 0.05;
X_des = [4,4]; %desired position
t_max = 60; %maximum simulation time
n = 4; %horizon
deltaT = 0.1;%timestep length
t = 1; %starting time
D_obs = 0.5;
alpha = 1;
h_values = zeros(1, t_max);
% v_fb_y_values = zeros(1, t_max);
% v_cmd_y_values = zeros(1,t_max);
while norm(X_des - X) > position_accuracy && t <= t_max
%describe behavior of obstacle------------------
% if t > 10 && t <= 25
% X_o(3) = -2;
% X_o(4) = -2;
% end
%
% if t > 25
% X_o(3) = 2;
% X_o(4) = 2;
% end
%-------------------------------------------------
v_fb = trajectoryController_h(X, X_des);
fprintf('TIME STEP: %d\n', t);
fprintf('v_fb: %d\n', v_fb);
fprintf('u value : %d\n', u);
[X, u] = safetyModule_h(X, u, v_fb, n, deltaT, X_o);
% v_fb_y_values(t) = v_fb(2);
% v_cmd_y_values(t) = u(2);
% %display h--------------------
d = sqrt((X(1) - X_o(1))^2 + (X(2) - X_o(2))^2);
% Calculate A and B
A = (X(1) - X_o(1)) / d;
B = (X(2) - X_o(2)) / d;
% Calculate constraint value for c(i)
h = - (A * u(1) + B * u(2) + alpha * (d - D_obs));
h_values(t) = h;
fprintf('h: %d \n \n ', h);
% %----------------------------
% fprintf('v_y: %d \n', u(2));
% %can use input to determine theta
theta = atan(u(2)/u(1));
% fprintf('Robot State: ');
% for i = 1:length(X)
% fprintf('%d ', X(i));
% end
% fprintf('\n \n');
%archive and plot
X_graph(t) = X(1);
Y_graph(t) = X(2);
cla;
daspect([1 1 1]);
hold on;
plot(X_o(1), X_o(2), 'or', 'MarkerSize', 5);
plot(X_des(1), X_des(2), 'og', 'MarkerSize', 5);
plot(X_graph(1:t), Y_graph(1:t), '-b'); %traveled path
% plot([X(1) X(1)+0.2*cos(theta)], [X(2) X(2)+0.2*sin(theta)], '-r'); %heading of robot
xlim([0, 4]);
ylim([0, 4]);
drawnow;
pause(deltaT);
%iterate
t = t + 1;
X_o(1) = X_o(1) + X_o(3)* deltaT;
X_o(2) = X_o(2) + X_o(4)* deltaT;
end
% figure;
% plot(1:t_max, h_values(1:t_max));
% xlabel('Timestep');
% ylabel('h Value');
% title('Value of h at Each Timestep');
%
% figure;
% plot(1:t_max, v_fb_y_values(1:t_max));
% hold on;
% plot(1:t_max, v_cmd_y_values(1:t_max));
% hold on;
% plot(1:t_max, v_cmd_y_values(1:t_max) - v_fb_y_values(1:t_max));
% xlabel('Timestep');
% ylabel('v');
% title('Value of obj at Each Timestep (only for v_y)');
% legend('v_{fb_y}', 'v_{cmd_y}', 'v_{cmd_y}-v_{fb_y}');
%{
safetyModule takes in system state, instantaneous input, feedback velocity, no. of
timesteps, timestep duration, no. of timesteps to return, and nearest obstacle's state
%}
% n discetized time points = N timesteps + 1
function [pos, u] = safetyModule_h(X, u, v_fb, n, deltaT, X_h_o)
%extend system state to match number of discretized time points in
%horizon including t = 0
X_h = [X, u, ones(1, 4 * (n-1))];
%initial guess (may not be feasible)
x0 = X_h;
%define constants
D_obs = 0.5; %margin of safety
alpha = 1;
%decompose v_fb
xdot_fb = v_fb(1);
ydot_fb = v_fb(2);
%objective function
fun = @(X_h) 0;
for i = 1:n
fun1 = @(X_h) abs(X_h(i+2*n) - xdot_fb) + abs(X_h(i+3*n) - ydot_fb);
fun = @(X_h) fun(X_h) + fun1(X_h);
end
%equality constraints
Aeq = zeros(2 * n + 2, 4 * n);
Aeq(1,1) = 1;
Aeq(2,1+n) = 1;
for i = 1:n-1
Aeq(2+i, i) = -1; Aeq(2+i,i+1) = 1; Aeq(2+i,i+n*2) = -deltaT;% x dynamics
Aeq(1 + n + i, n + i) = -1; Aeq(1 + n + i, n + i + 1) = 1; Aeq(1 + n + i, 3 * n + i) = -deltaT;% y dynamics
end
%initial velocity constraints
Aeq(2*n + 1, 2*n + 1) = 1;
Aeq(2*n + 2, 3*n+1) = 1;
beq = zeros(2*n + 2, 1);
beq(1)=X_h(1);
beq(2)=X_h(2);
beq(2*n + 1)=X_h(3);
beq(2*n + 2)=X_h(4);
% fprintf('Aeq:\n');
% for i = 1:size(Aeq, 1)
% for j = 1:size(Aeq, 2)
% fprintf('%d\t', Aeq(i, j)); % Change '%d' to '%f' for floating-point numbers
% end
% fprintf('\n'); % New line for each row
% end
%inequality constraints
nonlcon = @(X_h) CBFcon_h(X_h, n, deltaT, X_h_o, alpha, D_obs);
A = zeros(4*(n-1),4*n);
for i = 1:n-1
%less than or equal to acceleration constraint
A(i, 2*n + i) = -1; A(i, 2*n + i + 1) = 1;
A((n-1) + i, 3*n + i) = -1; A((n-1) + i, 3*n + i + 1) = 1;
%greater than or equal to acceleration
A(2*(n-1) + i, 2*n + i) = 1; A(2*(n-1) + i, 2*n + i + 1) = -1;
A(3*(n-1) + i, 3*n + i) = 1; A(3*(n-1) + i, 3*n + i + 1) = -1;
end
b = deltaT * 3 * ones(4*(n-1), 1);
% fprintf('beq:\n');
% for i = 1:size(beq, 1)
% for j = 1:size(beq, 2)
% fprintf('%d\t', beq(i, j)); % Change '%d' to '%f' for floating-point numbers
% end
% fprintf('\n'); % New line for each row
% end
%
% fprintf('Aeq:\n');
% for i = 1:size(Aeq, 1)
% for j = 1:size(Aeq, 2)
% fprintf('%d\t', Aeq(i, j)); % Change '%d' to '%f' for floating-point numbers
% end
% fprintf('\n'); % New line for each row
% end
%bounds
ub = [4 * ones(2*n, 1), 1.5 * ones(2*n, 1)];
lb = [zeros(2*n, 1), -1.5 * ones(2*n, 1)];
%run the solver
options = optimoptions('fmincon', 'Algorithm','sqp','Display', 'off','MaxIter',10000,'MaxFunEvals',10000);
[dec,~,exitflag] = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon, options);
exitflag
any(A*dec.'-b>1e-2)
any(abs(Aeq*dec.'-beq)>1e-2)
%Table
fields = {'t (horizon)', 'x', 'xdot', 'y', 'ydot'};
arrayCell = cell(n, numel(fields));
for i = 1:n
newArray = {i, dec(i), dec(i + 2 * n), dec(i + n), dec(i + 3 * n)};
arrayCell(i, :) = newArray;
end
T = cell2table(arrayCell, 'VariableNames', fields);
disp(T);
%extract x, y for the upcoming timestep
pos = zeros(1, 2);
pos(1) = dec(2);
pos(2) = dec(2+n);
%extract xdot, ydot for upcoming timestep
u = zeros(1, 2);
u(1) = dec(2*n + 2);
u(2) = dec(3*n + 2);
end
%{
trackingController calculates the desired velocity based on current,
position, desired position, and previous timestep input
%}
function v_fb = trajectoryController_h(X, X_des)
%define proportional constants for P controller
K_x = 1;
K_y = 1;
x = X(1);
y = X(2);
x_des = X_des(1);
y_des = X_des(2);
%extract errors
ex = x_des - x;
ey = y_des - y;
%Proportional controller
vfb_x = K_x * ex;
vfb_y = K_y * ey;
v_fb = [vfb_x vfb_y];
% final stretch constant speed definition:
% if (norm(X_des - x,y) < 0.75)
% mag = norm(v_fb);
% v_fb = (v_fb/mag) * 0.5;
% end
end
%{
function takes in robot system state, number of discretized time points in
the horizon, timestep duration, state of the nearest obstacle, and the CBF
constant
%}
%solver makes key assumption that obstacle does not accelerate
function [c, ceq] = CBFcon_h(X_h, n, deltaT, X_o, alpha, D_obs)
x_curr = X_o(1);
y_curr = X_o(2);
x_vel = X_o(3);
y_vel = X_o(4);
c = zeros(1, n);
% Loop through iterations
for i = 1:n
% Calculate distance between X_h(i) and current position
d = sqrt((X_h(i) - x_curr)^2 + (X_h(i + n) - y_curr)^2);
% Calculate A and B
A = (X_h(i) - x_curr) / d;
B = (X_h(i + n) - y_curr) / d;
% Calculate constraint value for c(i)
c(i) = - (A * X_h(i + 2*n) + B * X_h(i + 3*n) + alpha * (d - D_obs));
% Update current positions using velocity and time step
x_curr = x_curr + x_vel * deltaT;
y_curr = y_curr + y_vel * deltaT;
end
% Return empty equality constraint array
ceq = [];
end
Dhruv
on 27 Mar 2024
I changed
t_max = 25; %maximum simulation time
to
t_max = 60; %maximum simulation time
and
options = optimoptions('fmincon', 'Display', 'off');
dec = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon, options);
to
options = optimoptions('fmincon', 'Algorithm','sqp','Display', 'off','MaxIter',10000,'MaxFunEvals',10000);
[dec,~,exitflag] = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon, options);
exitflag
any(A*dec.'-b>1e-2)
any(abs(Aeq*dec.'-beq)>1e-2)
I don't understand your code - so I don't know what's going wrong when you set n = 10.
I personally find it difficult to comprehend code without comprehensive annotations explaining how the code functions and without seeing the mathematical principles underlying the algorithm's construction.
Based on my understanding, it appears that the program aims to calculate the steering angle for the particle to navigate from the starting point to the destination while avoiding obstacles, using a custom model predictive control (MPC) method inspired by artificial potential fields. Unfortunately, the code exhibits flaws when an obstacle is positioned in close proximity to the destination.
X = [0, 0]; % starting state
u = [0, 0]; % input at t = 0
% starting state of obstacle:
X_o = zeros(1,4);
X_o(1) = 3.5;
X_o(2) = 3.5;
X_des = [4, 4]; % desired position
Dhruv
on 27 Mar 2024
Dhruv
on 28 Mar 2024
Sam Chak
on 28 Mar 2024
It seems that your code fails to guide the particle to its destination when the destination is positioned within the repulsion radius of an obstacle. Troubleshooting complex codes can be challenging without detailed annotations that clarify the code's functionality and without insight into the underlying mathematical principles governing the algorithm's design.
If you're unable to identify the cause of one or more violated constraints, I recommend revisiting the fundamental mathematical equations and inequalities that form the basis of these constraints. It's crucial to examine the mathematics first and then review the code. Please ensure that you include detailed "comments" for each line in the constraint section, obstacle avoidance section, and steering section.
Here are the meanings of the exitflags from the documentation available under
exitflag — Reason fmincon stopped
integer
Reason fmincon stopped, returned as an integer.
All Algorithms:
1
First-order optimality measure was less than options.OptimalityTolerance, and maximum constraint violation was less than options.ConstraintTolerance.
0
Number of iterations exceeded options.MaxIterations or number of function evaluations exceeded options.MaxFunctionEvaluations.
-1
Stopped by an output function or plot function.
-2
No feasible point was found.
All algorithms except active-set:
2
Change in x was less than options.StepTolerance and maximum constraint violation was less than options.ConstraintTolerance.
trust-region-reflective algorithm only:
3
Change in the objective function value was less than options.FunctionTolerance and maximum constraint violation was less than options.ConstraintTolerance.
active-set algorithm only:
4
Magnitude of the search direction was less than 2*options.StepTolerance and maximum constraint violation was less than options.ConstraintTolerance.
5
Magnitude of directional derivative in search direction was less than 2*options.OptimalityTolerance and maximum constraint violation was less than options.ConstraintTolerance.
interior-point, sqp-legacy, and sqp algorithms:
-3
Objective function at current iteration went below options.ObjectiveLimit and maximum constraint violation was less than options.ConstraintTolerance.
Categories
Find more on Numerical Integration and Differentiation in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!




