Skip to content

Commit

Permalink
working pds, though differential_lift_x is _quite_ bad
Browse files Browse the repository at this point in the history
  • Loading branch information
apozharski committed May 24, 2024
1 parent 11f9305 commit c6b8d2e
Show file tree
Hide file tree
Showing 8 changed files with 227 additions and 19 deletions.
8 changes: 4 additions & 4 deletions +nosnoc/+dcs/Gcs.m
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ function generate_variables(obj, opts)
dims = obj.dims;
% dimensions
dims.n_lambda = dims.n_c;
obj.lambda = define_casadi_symbolic(opts.casadi_symbolic_mode,'lambda',obj.dims.n_lambda);
obj.lambda = define_casadi_symbolic(opts.casadi_symbolic_mode,'lambda',dims.n_lambda);

obj.dims = dims;
end
Expand All @@ -35,14 +35,14 @@ function generate_equations(obj, opts)

nabla_c = model.c.jacobian(model.x)';

obj.f_x = model.f_x + obj.E*nabla_c*obj.lambda;
obj.f_x = model.f_x + model.E*nabla_c*obj.lambda;

obj.f_x_fun = Function('f_x', {model.x, model.z, obj.lambda, model.u, model.v_global, model.p}, {obj.f_x, model.f_q});
obj.f_q_fun = Function('f_q', {model.x, model.z, obj.lambda, model.u, model.v_global, model.p}, {model.f_q});
obj.g_z_fun = Function('g_z', {model.x, model.z, model.u, model.v_global, model.p}, {model.g_z});
obj.g_alg_fun = Function('g_alg', {model.x, model.z, obj.lambda, model.u, model.v_global, model.p}, {[]});
obj.c_fun = Function('c_fun', {model.x, model.z, model.v_global, model.p}, {model.c})
obj.nabla_c_fun = Function('c_fun', {model.x, model.z, model.v_global, model.p}, {nabla_c})
obj.c_fun = Function('c_fun', {model.x, model.z, model.v_global, model.p}, {model.c});
obj.nabla_c_fun = Function('c_fun', {model.x, model.z, model.v_global, model.p}, {nabla_c});
obj.g_path_fun = Function('g_path', {model.x, model.z, model.u, model.v_global, model.p}, {model.g_path}); % TODO(@anton) do dependence checking for spliting the path constriants
obj.g_comp_path_fun = Function('g_comp_path', {model.x, model.z, model.u, model.v_global, model.p}, {model.g_comp_path});
obj.g_terminal_fun = Function('g_terminal', {model.x, model.z, model.v_global, model.p_global}, {model.g_terminal});
Expand Down
1 change: 0 additions & 1 deletion +nosnoc/+discrete_time_problem/Gcs.m
Original file line number Diff line number Diff line change
Expand Up @@ -304,7 +304,6 @@ function generate_direct_transcription_constraints(obj)
x_ijk = obj.w.x(ii,jj,opts.n_s+1);
z_ijk = obj.w.z(ii,jj,opts.n_s+1);
lambda_ijk = obj.w.lambda(ii,jj,opts.n_s+1);
mu_ijk = obj.w.mu(ii,jj,opts.n_s+1);

obj.g.dynamics(ii,jj,opts.n_s+1) = {x_ijk - x_ij_end};
obj.g.c_lb(ii,jj,opts.n_s+1) = {dcs.c_fun(x_ijk, z_ijk, v_global, p), 0, inf};
Expand Down
12 changes: 6 additions & 6 deletions +nosnoc/+model/Pds.m
Original file line number Diff line number Diff line change
Expand Up @@ -26,15 +26,15 @@ function verify_and_backfill(obj, opts)

dims.n_c = size(obj.c,1);

if ~isempty(E)
if size(E, 1) ~= size(E,2)
error("nosnoc: Projection matrix E must be square.");
if ~isempty(obj.E)
if size(obj.E, 1) ~= dims.n_x
error("nosnoc: Projection matrix E must be an n_x by n_x matrix where n_x is the number of functions defining the set C.")
end
if size(E, 1) ~= dims.n_c
error("nosnoc: Projection matrix E must be an n_c by n_c matrix where n_c is the number of functions defining the set C.")
if size(obj.E, 2) ~= dims.n_x
error("nosnoc: Projection matrix E must be an n_x by n_x matrix where n_x is the number of functions defining the set C.")
end
else
E = eye(dims.n_c);
obj.E = eye(dims.n_c);
end

obj.dims = dims;
Expand Down
6 changes: 5 additions & 1 deletion +nosnoc/+ocp/Solver.m
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,11 @@
case "nosnoc.model.Cls"
error("not implemented")
case "nosnoc.model.Pds"
error("not implemented")
obj.dcs = nosnoc.dcs.Gcs(model);
obj.dcs.generate_variables(opts);
obj.dcs.generate_equations(opts);
obj.discrete_time_problem = nosnoc.discrete_time_problem.Gcs(obj.dcs, opts);
obj.discrete_time_problem.populate_problem();
otherwise
error("Unknown model type")
end
Expand Down
27 changes: 20 additions & 7 deletions +nosnoc/Integrator.m
Original file line number Diff line number Diff line change
Expand Up @@ -39,17 +39,30 @@
obj.dcs.generate_equations(opts);
obj.discrete_time_problem = nosnoc.discrete_time_problem.Stewart(obj.dcs, opts);
obj.discrete_time_problem.populate_problem();
elseif opts.dcs_mode == DcsMode.Step % TODO: RENAME
error("not implemented")
elseif opts.dcs_mode == DcsMode.Step
obj.dcs = nosnoc.dcs.Heaviside(model);
obj.dcs.generate_variables(opts);
obj.dcs.generate_equations(opts);
obj.discrete_time_problem = nosnoc.discrete_time_problem.Heaviside(obj.dcs, opts);
obj.discrete_time_problem.populate_problem();
else
error("PSS models can only be reformulated using the Stewart or Heaviside Step reformulations.")
end
case "nosnoc.model.heaviside"
error("not implemented")
case "nosnoc.model.cls"
error("not implemented")
case "nosnoc.model.pds"
case "nosnoc.model.Heaviside"

obj.dcs = nosnoc.dcs.Heaviside(model);
obj.dcs.generate_variables(opts);
obj.dcs.generate_equations(opts);
obj.discrete_time_problem = nosnoc.discrete_time_problem.Heaviside(obj.dcs, opts);
obj.discrete_time_problem.populate_problem();
case "nosnoc.model.Cls"
error("not implemented")
case "nosnoc.model.Pds"
obj.dcs = nosnoc.dcs.Gcs(model);
obj.dcs.generate_variables(opts);
obj.dcs.generate_equations(opts);
obj.discrete_time_problem = nosnoc.discrete_time_problem.Gcs(obj.dcs, opts);
obj.discrete_time_problem.populate_problem();
end
end

Expand Down
47 changes: 47 additions & 0 deletions examples/new_design_examples/new_design_moving_set.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,47 @@

clear all
clc
close all
import casadi.*
%% discretization parameters
N_sim = 100;
T_sim = 10;

%% NOSNOC settings
problem_options = nosnoc.Options();
solver_options = nosnoc.solver.Options();
problem_options.n_s = 2;
settings.homotopy_update_slope = 0.1;
%problem_options.rk_scheme = RKSchemes.GAUSS_LEGENDRE;
problem_options.rk_scheme = RKSchemes.RADAU_IIA;
problem_options.rk_representation= 'differential_lift_x';
%problem_options.rk_representation= 'integral';
problem_options.N_sim = N_sim;
problem_options.N_finite_elements = 2;
problem_options.T_sim = T_sim;

%solver_options.homotopy_steering_strategy = 'ELL_INF';


model = nosnoc.model.Pds();

model.x0 = [0;pi-2;0];
x = SX.sym('x',2);
t = SX.sym('t');
model.x = [x;t];
model.c = [-norm(x - [sin(t);cos(t)])^2+(1-0.5)^2];
model.f_x = [0; 0; 1];
model.E = diag([1,1,0]);

integrator = nosnoc.Integrator(model, problem_options, solver_options);
[t_grid, x_res] = integrator.simulate();
%
figure
plot(x_res(1,:), x_res(2,:))
grid on
xlabel('$t$','Interpreter','latex')
ylabel('$x(t)$','Interpreter','latex')
grid on

fig = figure('Position', [10 10 1600 800]);
plot_moving_set([],x_res,[0.5], ["circle"], fig, 'moving_set');
44 changes: 44 additions & 0 deletions examples/new_design_examples/new_design_pds.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,44 @@

clear all
clc
close all
import casadi.*
%% discretization parameters
N_sim = 100;
T_sim = 10;

%% NOSNOC settings
problem_options = nosnoc.Options();
solver_options = nosnoc.solver.Options();
problem_options.n_s = 2;
settings.homotopy_update_slope = 0.1;
%problem_options.rk_scheme = RKSchemes.GAUSS_LEGENDRE;
problem_options.rk_scheme = RKSchemes.RADAU_IIA;
problem_options.rk_representation= 'differential_lift_x';
%problem_options.rk_representation= 'integral';
problem_options.N_sim = N_sim;
problem_options.N_finite_elements = 2;
problem_options.T_sim = T_sim;
solver_options.homotopy_steering_strategy = 'ELL_INF';


model = nosnoc.model.Pds();

model.x0 = [sqrt(2)/2;sqrt(2)/2];
x = SX.sym('x',2);
model.x = x;
model.c = x(2)+0.2;
model.f_x = [x(2);-x(1)];

model.x0 = [0;pi-2;0];

integrator = nosnoc.Integrator(model, problem_options, solver_options);
[t_grid, x_res] = integrator.simulate();
%
figure
plot(x_res(1,:), x_res(2,:))
grid on
xlabel('$t$','Interpreter','latex')
ylabel('$x(t)$','Interpreter','latex')
grid on

101 changes: 101 additions & 0 deletions examples/new_design_examples/plot_moving_set.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,101 @@
function plot_moving_set(h,pos,r,type,fig,vidname)
n = length(r);
if ~exist('type')
type = [];
for ii=1:n
type = [type, "circle"];
end
end
if length(type) ~= n
error('not all types specified')
end
indices = {};
curr_idx = 1;
for ii=1:n
switch type(ii)
case "circle"
indices{ii} = curr_idx:curr_idx+1;
curr_idx = curr_idx+2;
case "square"
indices{ii} = curr_idx:curr_idx+2;
curr_idx = curr_idx+3;
end
end
if ~exist('fig')
fig = figure;
else
figure(fig)
end
axis equal
box on;
ax = gca;
set(ax,'XTick',[-10,0,10], 'YTick', [0,5,10], 'TickLength', [0.025,0.025], 'LineWidth', 4.0);
xlim([-3,3])
ylim([-3,3])
ax.XAxis.FontSize = 36;
ax.YAxis.FontSize = 36;
xlabel("$c_x$", "fontsize", 32)
ylabel("$c_y$", "fontsize", 32)

%axis off
discs = {};
tau = linspace(0, 2*pi)';
rot = @(theta) [cos(theta) -sin(theta);...
sin(theta) cos(theta)];
circ = @(p,r) [r*cos(tau)+p(1),r*sin(tau)+p(2)];
sqr = @(p,r) transpose(rot(-p(3))*([r,r;r,-r;-r,-r;-r,r]') + p(1:2));
for ii=1:n
p = pos(indices{ii},1);
switch type(ii)
case "circle"
v = circ(p,r(ii));
case "square"
v = sqr(p,r(ii));
end
if ii==3
color = [0.8500 0.3250 0.0980];
else
color = [0 0.4470 0.7410];
end
discs{ii} = patch(v(:,1),v(:,2), color);
end
t = pos(end,1);
v_C = circ([sin(t);cos(t)],1);
C = patch(v_C(:,1),v_C(:,2), [1 0 0], 'FaceAlpha', 0.0, 'LineStyle', '--', 'EdgeColor' , [1 0 0]);
if exist('vidname')
mkdir([vidname '_frames'])
writer = VideoWriter([vidname '.avi']);
open(writer);
frame = getframe(gca);
exportgraphics(gca, [vidname '_frames/' num2str(1) '.pdf'])
writeVideo(writer,frame)
end
pause(0.1);
for jj=2:(size(pos,2)-1)
for ii=1:n
p = pos(indices{ii},jj+1);
switch type(ii)
case "circle"
v = circ(p,r(ii));
case "square"
v = sqr(p,r(ii));
end
discs{ii}.Vertices = v;
end
t = pos(end,jj);
v_C = circ([sin(t);cos(t)],1);
C.Vertices = v_C;
drawnow limitrate;
%pause(h(jj));
pause(0.1)
if exist('vidname')
frame = getframe(gca);
ff=jj+1;
%exportgraphics(gca, [vidname '_frames/' num2str(jj) '.pdf'])
writeVideo(writer,frame);
end
end
if exist('vidname')
close(writer);
end
end

0 comments on commit c6b8d2e

Please sign in to comment.