The boundary value problem is

$$ \begin{cases} \dot{x}_i = \begin{cases} (0.5D^{-1}\psi)_i, \text{ if }(0.5D^{-1}\psi)_i \le 0 \\ 0 \text{, otherwise} \end{cases} \\ \dot{\psi} = 2\Sigma x \\ x(0) = x^0 \\ x(T) = 0 \end{cases} $$

Here $x^0>0$, $T$ is given, $D>0$ is a diagonal matrix, $\Sigma>0$ is a symmetric matrix. Actually the constraint is $x(t)\ge 0\forall t$, equivalent to $\dot{x}(t) \le 0$.

For a 2-dimension case, I try to solve it in MATLAB with this code

function [tint, traj] = liq2mon_v3(x0, D, sigma, T)

  tlow = 0;
  thigh = T;
  L = length(x0); % L=2 in this example
  nsegm = 100; % I know it is for final calc only, not for the mesh
  solinit = bvpinit(linspace(tlow, thigh, nsegm), @guess);
  sol = bvp5c(@bvp4ode, @bvp4bc, solinit);
  tint = linspace(tlow, thigh);
  z = deval(sol, tint);
  traj = z(1:L, :);

  function xpsi = guess(t)
    xpsi = zeros(1, L*2);
    xpsi(1:L) = (x0*(1.0 - 1.0/T*t))';

  function dydt = bvp4ode(t, y)
    xcomp = y(1:L); % L=2 in this example
    psicomp = y(L+1:end);
    u    = D \ psicomp / 2.0;
    dx   = min(u, 0); 
    dpsi = 2.0 * sigma * xcomp;
    dydt = [dx' dpsi']; 

  function res = bvp4bc(ya, yb)
    xcompa = ya(1:2);
    xcompb = yb(1:2);
    res = [(xcompa - x0)' (xcompb)'];


and I get an error

*Unable to solve the collocation equations -- a singular Jacobian  encountered
Error in ==> liq2mon_v3 at 13
sol = bvp5c(@bvp4ode, @bvp4bc, solinit)*

However if I change the line

dx = min(u, 0);


u(u>0) = 1e-13;
dx = u;

I do get a solution, with $\min(\min(x)) \approx -8\times 10^{-14}$ for the same $x^0$.


  1. Is there a way to solve this system in MATLAB (or another app) so that the constraint is not violated?
  2. Maybe this could be solved in Python, even with the current constraint violation?
  3. Why is the Jacobian singular somewhere within MATLAB solver in this case?
  • $\begingroup$ I believe now I should add a comment, if not an 'answer'. This problem comes from LQR (MPC) with $F=x'\Sigma x + \dot{x}Dx$. As I learned recently from the literature (for instance, "Tutorial Overview of Model Predictive Control"), the target state should be within the allowed region for both control and state whereas in this problem it was on the boundary. 'Relaxing' the constraint to $10^{-13}$ made the origin the internal point and everything worked according to the books. However it is still unclear what exactly went wrong in MATLAB solver. $\endgroup$ Oct 30 '15 at 10:08
  • $\begingroup$ Of course $F=x'\Sigma x + \dot{x'}D\dot{x}$. I ran out of time to edit the first comment. $\endgroup$ Oct 30 '15 at 10:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.