Difference between revisions of "Cushioned Oscillation (PROPT)"
From mintOC
FelixMueller (Talk | contribs) |
FelixMueller (Talk | contribs) |
||
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]] to create the reference solution and its plot for the [[Cushioned Oscillation]] problem. |
<source lang="matlab"> | <source lang="matlab"> |
Revision as of 18:49, 31 January 2016
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.
%% 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)