[n,m] = size(J_ux);
  if (n != n_U)||(m != n_x)
    error("J_ux should be %ix%i not %ix%i",n_U,n_x,n,m);
  endif


  if length(gamma)>0		# Constraints exist: do the QP algorithm
    ## QP solution for weights U	
    [U,iterations] = qp_mu(J_uu,(J_ux*x - J_uw*W),Gamma,gamma,mu); # QP solution for weights U	
    [U,iterations] = qp_mu(J_uu,(J_ux*x - J_uw*W),Gamma,gamma,mu);

    ##U = qp(J_uu,(J_ux*x - J_uw*W),Gamma,gamma); # QP solution for weights U
    ##U = pd_lcp04(J_uu,(J_ux*x - J_uw*W),Gamma,gamma); # QP solution for weights U
    u = Us0*U;			# Control signal
  else			# Do the unconstrained solution
    ## Compute the open-loop gains
    iterations = 0;