Difference between revisions of "Cushioned Oscillation (PROPT)"

From mintOC
Jump to: navigation, search
 
(2 intermediate revisions by the same user not shown)
Line 1: Line 1:
Below you can find the MATLAB file that was used in [[Category: Tomdyn/PROPT]] to create the reference solution and its plot for the [[Cushioned Oscillation]] problem.
+
Below you can find the MATLAB file that was used in [[:Category: TomDyn/PROPT | TomDyn/PROPT]] to create the reference solution and its plot for the [[Cushioned Oscillation]] problem.
  
 
<source lang="matlab">
 
<source lang="matlab">

Latest revision as of 18:49, 31 January 2016

Below you can find the MATLAB file that was used in TomDyn/PROPT to create the reference solution and its plot for the Cushioned Oscillation problem.

%% Cushioned Oscillation
% (c) Maximilian von Wolff
% 
%% Problem Setup
%
%
%
%Problem Parameters
t_0  = 0;
 
x_0 = 2;  %starting position
v_0 = 5;  %starting velocity in m/s
umm = 5;  %control constraint
m   = 5; %mass in kg
c   = 10;  %spring stiffness in N/m
 
 
n=80; %Number of collocation points
 
 
%Setup
 
toms t t_f            %independent variables
p = tomPhase('p', t, t_0, t_f, n);
setPhase(p);
tomStates x v        
tomControls u
 
%initial states
xi = [x_0; v_0];
 
%initial guess
x0 = {t_f == 10
    collocate({u == 0
                 })
     icollocate({x == xi(1)
                 v == xi(2)
                 })};
 
%Box constraints
cbox = {-umm <= collocate(u) <= umm
        };
 
%Boundary constraints
cbnd = {initial({   x == xi(1); 
                    v == xi(2);
                    })
        final({     x == 0;
                    v == 0;
                    })};
 
% ODE's
 
dx  = v;
dv  = 1./m.*(u-c*x);
 
 
ceq = collocate({
        dot(x) == dx
        dot(v) == dv
        });
 
 
objective = t_f;    
 
%% Solve the problem
    options = struct;
    options.name = 'cushioned oscillation';
    solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);
 
    tfopt = subs(t_f,solution);
    xopt = subs(x,solution);
    vopt = subs(v,solution);
    uopt = subs(u,solution);
 
 
%Plotting solution
 
figure(1)
subplot(3,1,1);
ezplot(x); legend('x');
xlabel('time');
ylabel('position in m');
title('Position');
 
subplot(3,1,2);
ezplot(v); legend('v');
xlabel('time');
ylabel('velocity in m/s');
title('Velocity');
 
subplot(3,1,3);
ezplot(u); legend('u');
xlabel('time');
title('Control');
 
disp('Final Time');
disp(solution.t_f)