Cushioned Oscillation (PROPT)

Below you can find the MATLAB file that was used in 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)