Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
fmincon not working for delta V calculations

Subject: fmincon not working for delta V calculations

From: Shen

Date: 30 Apr, 2009 16:49:01

Message: 1 of 1

I'm using MATLAB 7.5.0 (R2007b). I'm trying to use fmincon to do a simple project. What I have is a planar earth to moon trajectory from a circular LEO to a circular LMO. There's only two delta V changes with one for LEO departure and the other for LMO insertion. I'm trying to minmize delta V assuming the time of flight is a constant 3 days. My state equations are the x, y positions and the u,v velocities. Unfortunately, when I run the code, only 1 iteration is done and the exitflag is -2 which means it failed to converge. Can someone tell me whether this is just a programming issue or whether I used a wrong constraint or is missing information? I set delVLEO and delVLMO to be global variables. Could that be a reason?

% Main Function
clear all
clc
global N; % Number of timesteps
global h; % step size
global delVLEO; global delVLMO;
N;
t0=0; tf$.*3600.*3; % Final time set to 3 days
h=(tf-t0)./N;
delVLEO=2; % Guess
delVLMO=2; % Guess

% Initial Guess
% Inertial reference frame contained in the Moon orbital plane.
% Origin is Earth center.
% x-axis points toward the Moon initial position
% y-axis perpendicular to x and pointed towards Moon's initial inertial v
VLEO = 7.633; % km/s
for i=1:N+1,
    xP(i) = 6841; % x position
    yP(i) = 0; % y position
    uP(i) = 0; % x velocity
    vP(i) = VLEO; % y velocity
end
% delVLEO = 0;
% delVLMO = 0;

x0(1:N+1) = xP; % This definition is just so that fmincon works.
x0((N+1)+1:2.*(N+1)) = yP;
x0(2.*(N+1)+1:3.*(N+1)) = uP;
x0(3.*(N+1)+1:4.*(N+1)) = vP;
% x0(4.*(N+1)+1) = delVLEO;
% x0(4.*(N+1)+2) = delVLMO;

% State Constraints
xPlb = -inf; xPub = +inf;
yPlb = -inf; yPub = +inf;
uPlb = -inf; uPub = +inf;
vPlb = -inf; vPub = +inf;
% delVLEOlb = 0; delVLEOub = 7; % km/s
% delVLMOlb = -7; delVLMOub = 0; % km/s

lb(1:N+1) = xPlb; % Definition just so fmincon works
lb((N+1)+1:2.*(N+1)) = yPlb;
lb(2.*(N+1)+1:3.*(N+1)) = uPlb;
lb(3.*(N+1)+1:4.*(N+1)) = vPlb;
% lb(4.*(N+1)+1) = delVLEOlb;
% lb(4.*(N+1)+2) = delVLMOlb;

ub(1:N+1) = xPub; % Definition just so fmincon works
ub((N+1)+1:2.*(N+1)) = yPub;
ub(2.*(N+1)+1:3.*(N+1)) = uPub;
ub(3.*(N+1)+1:4.*(N+1)) = vPub;
% ub(4.*(N+1)+1) = delVLEOub;
% ub(4.*(N+1)+2) = delVLMOub;

% No equality constraints
A = [];
b = [];
Aeq = [];
beq = [];

% Define the constraints function to use
funcon = 'funcon_rungekutta';
% Define the cost function to use
funobj = 'funobj_trap';

options=[];
options=optimset(options,'GradObj','off','GradConstr','off','display','off','LargeScale','off');
tic
[x, f, inform, output, lambda, g, H] = fmincon(funobj, x0, A, b, Aeq, beq, lb, ub, funcon, options)
toc



%! CONSTRAINT FUNCTION
function [cin,c]=funcon_trap(x)
global N;
global h;
global delVLEO; global delVLMO;
% display('hi!')
% Define constants
muE = 3.986E5; % km3/s2 earth gravitational constant
muM = 4.903E3; % km3/s2 lunar grav. constant
rME = 3.844E5; % km moon radial distance from Earth center
omegaM = 2.6491E-6; % rad/s moon angular velocity
VM = 1.0183; % km/s moon inertial velocity
RE = 6378; % km radius of earth
RM = 1738; % km radius of moon

rLEO = 6841; % km
VLEO = 7.633; % km/s
rLMO = 1838; % km
VLMO = 1.633; % km/s

% Do some reassignment
xP = x(1:N+1);
yP = x((N+1)+1:2.*(N+1));
uP = x(2.*(N+1)+1:3.*(N+1));
vP = x(3.*(N+1)+1:4.*(N+1));
% delVLEO = x(4.*(N+1)+1);
% delVLMO = x(4.*(N+1)+2);
% display('hi2!')
%! Define more variables
rPE = sqrt(xP.^2 + yP.^2); % radial distance of s/c from Earth
for i=1:N+1,
    t=h.*(i-1);
    thetaM(i) = omegaM.*t; % angular coordinate associated with the moon position
end

xM = rME.*cos(thetaM); % moon x-position
yM = rME.*sin(thetaM); % moon y-position
rPM = sqrt((xP-xM).^2+(yP-yM).^2);% radial distance of s/c from moon
uM = -VM.*sin(thetaM); % moon orbit velocity along x
vM = VM.*cos(thetaM); % moon orbit velocity along y
% relative-to-moon cordinates
xPM = xP - xM;
yPM = yP - yM;
uPM = uP - uM;
vPM = vP - vM;

% Write out the definitions!
Xdot(1,:) = uP; % xPdot 1st row is just a definition
Xdot(2,:) = vP; % yPdot 2nd row also
Xdot(3,:) = -(muE./rPE.^3).*xP - (muM./rPM.^3).*(xP-xM); % uPdot 3rd row is orbital mechanics
Xdot(4,:) = -(muE./rPE.^3).*yP - (muM./rPM.^3).*(yP-yM); % vPdot 4th row is orbital mechanics

% Equality constraints. (Just the state equations!!)
% Just using a very rough step increment approximation
c(1:N) = xP(1:end-1)-xP(2:end) + h.*uP(1:end-1);
c(N+1:2.*N) = yP(1:end-1)-yP(2:end) + h.*vP(1:end-1);
c(2.*N+1:3.*N) = uP(1:end-1)-uP(2:end) + h.*Xdot(3,1:end-1);
c(3.*N+1:4.*N) = vP(1:end-1)-vP(2:end) + h.*Xdot(4,1:end-1);

% Boundary conditions (departure conditions. velocity direction, radius,
% velocity magnitude)
c(4.*N+1)=xP(1).*uP(1)+yP(1).*vP(1); % This must be 0 for circular orbit
c(4.*N+2)=xP(2).*uP(2)+yP(2).*vP(2); % This must be 0 for accelearting velocity impulse to be tangential to LEO
c(4.*N+3)=rLEO-sqrt(xP(1).^2 + yP(1).^2); % Initial radius is at circular LEO.
c(4.*N+4)=VLEO-sqrt(uP(1).^2 + vP(1).^2); % Initial velocity is just circular tangential velocity at LEO
% More boundary conditions. Arrival conditions.
c(4.*N+5)=rLMO-sqrt((xP(end-1)-xM(end-1)).^2+(yP(end-1)-yM(end-1)).^2); % arrival radius must be at low moon orbit
c(4.*N+6)=rLMO-sqrt((xP(end)-xM(end)).^2+(yP(end)-yM(end)).^2); % final radius must still be at LMO
c(4.*N+7)=VLMO-sqrt((uP(end)-uM(end)).^2+(vP(end)-vM(end)).^2); % final velocity is just circular tang. vel. at LMO
c(4.*N+8)=xPM(end-1).*uPM(end-1)+yPM(end-1).*vPM(end-1); % tangential requirements
c(4.*N+9)=xPM(end).*uPM(end)+yPM(end).*vPM(end); % tangential requirements

% Inequality constraints (must be less than 0 obviously)
cin = [];
cin(1) = -delVLEO;
cin(2) = -delVLMO; % deltaVs must always be positive




%! COST FUNCTION
function [objf] = funobj_trap(x)
global N;
global h;
global delVLEO; global delVLMO;
% relative-to-earth coordinates
% xP = x(1:N+1);
% yP = x((N+1)+1:2.*(N+1));
uP = x(2.*(N+1)+1:3.*(N+1));
vP = x(3.*(N+1)+1:4.*(N+1));
% moon motion
% rME = 3.844E5; % km moon radial distance from Earth center
omegaM = 2.6491E-6; % rad/s moon angular velocity
VM = 1.0183; % km/s moon inertial velocity
for i=1:N+1,
    t=h.*(i-1);
    thetaM(i) = omegaM.*t; % angular coordinate associated with the moon position
end
% xM = rME.*cos(thetaM); % moon x-position
% yM = rME.*sin(thetaM); % moon y-position
% rPM = sqrt((xP-xM).^2+(yP-yM).^2);% radial distance of s/c from moon
uM = -VM.*sin(thetaM); % moon orbit velocity along x
vM = VM.*cos(thetaM); % moon orbit velocity along y
% relative-to-moon coordinates
% xPM = xP - xM;
% yPM = yP - yM;
uPM = uP - uM;
vPM = vP - vM;

delVLEO = sqrt(uP(2).^2+vP(2).^2)-sqrt(uP(1).^2+vP(1).^2); % tangential so just subtract magnitudes
delVLMO = sqrt(uPM(end).^2+vPM(end).^2)-sqrt(uPM(end-1).^2+vPM(end-1).^2);

% Define cost function
objf = delVLEO+delVLMO;

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us