Cushioned Oscillation (PROPT)

Below you can find the MATLAB file that was used to create the reference solution and its plot.

%% 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)