search

casadi / casadi

CasADi is a symbolic framework for numeric optimization implementing automatic differentiation in forward and reverse modes on sparse matrix-valued computational graphs. It supports self-contained C-code generation and interfaces state-of-the-art codes such as SUNDIALS, IPOPT etc. It can be used from C++, Python or Matlab/Octave.
http://casadi.org
GNU Lesser General Public License v3.0
1.71k stars 379 forks source link

Bug in mex interface #2375

Closed jgillis closed 5 years ago

jgillis commented 5 years ago 
import casadi.*

T = 10; % Time horizon
N = 2; % number of control intervals

F = @(x,u) x;
opti = casadi.Opti();

x = opti.variable(2,N+1);
u = opti.variable(1,N);

opti.minimize(sumsqr(x));

for k=1:N
  opti.subject_to(x(:,k+1)==F(x(:,k),u(:,k)));
end
opti.subject_to(-1<=u<=1);

p = opti.parameter(2);
opti.subject_to(x(:,1)==p);

opti.solver('sqpmethod',struct('qpsol','qrqp'));

opti.set_value(p,[0;1]);
sol = opti.solve();
res = to_function(opti);
x = res.x;
M = casadi.Function('M',{p},{x(2*(N+1)+1)});

M.generate('Mc',struct('mex',true,'main',true));
mex Mc.c -largeArrayDims

Mc('M',[-1;-1])
function res = to_function(opti)
import casadi.*
sp = sparse(DM(sparsity(jacobian(opti.g,opti.x)),1));

is_single = full(sum(sp,2)==1);
is_linear = ~which_depends(opti.g,opti.x,2,true)';
is_simple = is_single & is_linear;

is_simple = is_simple & false;
[col,~] = find(sp(is_simple,:)');
%assert(numel(unique(col))==numel(col))

g = opti.g;

% Detect  f2(p)x+f1(p)==0
gf = Function('gf',{opti.x,opti.p},{g(find(is_simple)),jtimes(g(find(is_simple)),opti.x,ones(opti.nx,1))});

[f1,f2] = gf(0,opti.p);

lbg = opti.lbg;
ubg = opti.ubg;

lbx = -MX.inf(opti.nx,1);
ubx = MX.inf(opti.nx,1);

lb = (lbg(find(is_simple))-f1)./abs(f2);
ub = (ubg(find(is_simple))-f1)./abs(f2);

lbc = {};
ubc = {};
for i=1:numel(col)
   j = col(i);
   lbx(j) = max(lbx(j),lb(i));
   ubx(j) = min(ubx(j),ub(i));
end

lbx(col) = (lbg(find(is_simple))-f1)./abs(f2);
ubx(col) = (ubg(find(is_simple))-f1)./abs(f2);

%%

solver_s = nlpsol('solver','sqpmethod',struct('x',opti.x,'p',opti.p,'f',opti.f,'g',g(find(~is_simple))),struct('qpsol','qrqp','error_on_fail',true,'qpsol_options',struct('print_iter',false,'print_header',false),'print_time',false,'print_iteration',false));
%solver_s = nlpsol('solver','sqpmethod',struct('x',opti.x,'p',opti.p,'f',opti.f,'g',g(find(~is_simple))),struct('qpsol','qrqp'));

res = solver_s('p',opti.p,'lbx',lbx,'ubx',ubx,'lbg',lbg(find(~is_simple)),'ubg',ubg(find(~is_simple)));
end
jgillis commented 5 years ago

Probable cause: https://github.com/casadi/casadi/commit/d1d1601b6e1db8cf73dd0a184435f3b93de3da9a#diff-e5faf994e3946dbcf42e5d7661cf0877R2157