Understanding how to add constraints

[G Prathap]

I have formulated the problem in the following way, 

In order to solve with Mosek, I have rearranged in the following way, 

I believe my reformulation from (1) to (2)  is correct.  But I have a problem with adding the following constraints, i.e, how to subtract and get l1 norm 

Here is the code that I have so far. 

#include <iostream>

#include "fusion.h"

using namespace mosek::fusion;

using namespace monty;

int main(int argc, char ** argv)

{

int number_of_steps = 6;

int lamda1 = 1000;

int lamda2 = 1000;

int dim_space = 3;

auto A = new_array_ptr<double, 2>({ {-0.470233, 0.882543 , 0},

{0.470233, -0.882543 , -0},

{-0.807883, -0.430453 , 0.402535},

{0.807883 , 0.430453 ,-0.402535},

{-0.355254 ,-0.189285 ,-0.915405},

{0.355254 , 0.189285 , 0.915405}});

auto b = new_array_ptr<double, 2>({{ 3.21407},

{ 0.785929},

{4.81961},

{-1.40824},

{0.878975},

{1.12103}});

auto start_ = new_array_ptr<double, 2>({{-4}, {0}, {2}});

auto end_ = new_array_ptr<double, 2>({{-1}, {-1}, {-1}});

Model::t M = new Model();

auto x = M->variable(new_array_ptr<int,1>({dim_space, number_of_steps}), Domain::unbounded());

auto t = M->variable(number_of_steps, Domain::unbounded());

// std::vector<Expression::t> diff_list(number_of_steps);

for(int i=1; i<number_of_steps-1; i++){

Variable::t x_i = Var::vstack(x->index(0,i), x->index(1,i), x->index(2,i));

Variable::t x_i1 = Var::vstack(x->index(0,i+1), x->index(1,i+1), x->index(2,i+1));

// How to define the difference between x_i1 - x_i and get the l1 norm, i.e., |x_i1 - x_i|

}

}

Could you help me with how to add those constraints properly?

Thanks,

Geesara 

[G Prathap]

After spending some time, I got something but does not give me the expected result. Here I am adding the original MATLAB version and C++ version I am having at the moment, problem is bit different that I've defined above

  

   Matlab code 

    clear all
    close all
    clc

    number_of_steps = 15;
    lambda3_val = 1000;
    lambda1_val = 1000;
    lambda2_val = 0.1;
    dim_ = 3;
    Ab = [-0.470233  0.882543         0   3.21407
     0.470233 -0.882543        -0  0.785929
    -0.807883 -0.430453  0.402535   4.81961
     0.807883  0.430453 -0.402535  -1.40824
    -0.355254 -0.189285 -0.915405  0.878975
     0.355254  0.189285  0.915405   1.12103];

    A = Ab(:,1:dim_);
    b = Ab(:,dim_+1);

    start = [-4 , 0, 2];
    goal = [-1, -1, 1];

    nPolytopes = 1;
    free_space = polytope(A, b);

    cvx_solver mosek
    cvx_begin
        variable x(3, number_of_steps)
        binary variable c(nPolytopes, number_of_steps);
        cost = 0;
        for i = 1:(number_of_steps-1)
            cost = cost + pow_pos(norm(x(:, i) - x(:, i+1), 1),2)*lambda2_val;
        end
        cost = cost + pow_pos(norm(x(:, 1) - start'), 2)*lambda1_val;
        cost = cost + pow_pos(norm(x(:, number_of_steps) - goal'),2)*lambda3_val;
        minimize(cost)
        subject to
            for i = 1:number_of_steps
                A*x(:, i) <= b;
            end
    cvx_end

Here is the code the conversion of the above into c++ using Mosek

    #include <iostream>
    #include "fusion.h"

    using namespace mosek::fusion;
    using namespace monty;

    int main(int argc, char ** argv)
    {

      int number_of_steps = 15;
      int lambda1_val = 1000;
      int lambda2_val = 1000;
      int lambda3_val = 0.1;
      int dim_space = 3;
      auto A = new_array_ptr<double, 2>({{-0.4702,    0.8825,         0},
                                          {0.4702,   -0.8825,         0},
                                          {-0.8079,   -0.4305,    0.4025},
                                          {0.8079,    0.4305,   -0.4025},
                                          {-0.3553,   -0.1893,   -0.9154},
                                          {0.3553,    0.1893,    0.9154}});
     
      auto b = new_array_ptr<double, 1>({3.2141,
                                        0.7859,
                                        4.8196,
                                      -1.4082,
                                        0.8790,
                                        1.1210});

      auto end_ = new_array_ptr<double, 1>({-1, -1, -1});
      auto start_ = new_array_ptr<double, 1>({-4, 0, 2});

      Model::t M = new Model();
      auto x = M->variable(new_array_ptr<int,1>({dim_space, number_of_steps}), Domain::unbounded()) ;
      auto t = M->variable(number_of_steps, Domain::unbounded());

      auto t_tmp = M->variable(number_of_steps, Domain::unbounded());
      auto lambda_1 = M->parameter("lambda_1");
      auto lambda_2 = M->parameter("lambda_2");
      auto lambda_3 = M->parameter("lambda_3");
     
      Variable::t x_0 = Var::vstack(x->index(0, 0),  x->index(1, 0), x->index(2, 0));
      M->constraint(Expr::vstack(t_tmp->index(0), Expr::sub(x_0, start_)), Domain::inQCone());
      M->constraint(Expr::hstack(t->index(0), 1, t_tmp->index(0)), Domain::inPPowerCone(1.0/2));

      for(int i=1; i<number_of_steps-1; i++){
        Variable::t x_i = Var::vstack(x->index(0,i),  x->index(1,i), x->index(2,i));
        Variable::t x_i1 = Var::vstack(x->index(0,i+1),  x->index(1,i+1), x->index(2,i+1));

        M->constraint(Expr::vstack(t_tmp->index(i), Expr::sub(x_i1, x_i)), Domain::inQCone());
        M->constraint(Expr::hstack(t->index(i), 1, t_tmp->index(i)), Domain::inPPowerCone(1.0/2));
      }

      Variable::t x_n = Var::vstack(x->index(0, number_of_steps-1),  x->index(1, number_of_steps-1), x->index(2, number_of_steps-1));
      M->constraint(Expr::vstack(t_tmp->index(number_of_steps-1), Expr::sub(x_n, end_)), Domain::inQCone());
      M->constraint(Expr::hstack(t->index(number_of_steps-1), 1, t_tmp->index(number_of_steps-1)), Domain::inPPowerCone(1.0/2));

      for (int i = 0; i < number_of_steps; i++)
      {
        auto x_i = Var::vstack(x->index(0,i), x->index(1, i), x->index(2, i));
        M->constraint(Expr::mul(A, x_i), Domain::lessThan(b));
      }

      M->setLogHandler([](const std::string& msg){std::cout<< msg << std::flush;});
   
      auto lambda1 = M->getParameter("lambda_1");
      auto lambda2 = M->getParameter("lambda_2");
      auto lambda3 = M->getParameter("lambda_3");
      lambda1->setValue(lambda1_val);
      lambda2->setValue(lambda2_val);
      lambda3->setValue(lambda3_val);
                 
      auto objs = new_array_ptr<Expression::t, 1>(number_of_steps);
      (*objs)[0] = (Expr::mul(lambda1, t->index(0)));

      for(int i=1; i<number_of_steps-1; i++){

        (*objs)[i] = Expr::mul(lambda2, t->index(i));
      }
      (*objs)[number_of_steps-1] = Expr::mul(lambda3, t->index(number_of_steps-1));

      M->objective(ObjectiveSense::Minimize, Expr::add(objs));

      M->solve();
      auto sol = *(x->level());
      std::cout<< "solution "<< sol << std::endl;

    }

Yet I won't get the result that I get from Matlab implementation, I am not sure `pow_pos(norm(x(:, 1) - start'), 2)` I have converted correctly, here is the converted code

          Variable::t x_0 = Var::vstack(x->index(0, 0),  x->index(1, 0), x->index(2, 0));
          M->constraint(Expr::vstack(t_tmp->index(0), Expr::sub(x_0, start_)), Domain::inQCone());
          M->constraint(Expr::hstack(t->index(0), 1, t_tmp->index(0)), Domain::inPPowerCone(1.0/2));

 



 

--
You received this message because you are subscribed to a topic in the Google Groups "mosek" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/mosek/35JseJyY8cU/unsubscribe.
To unsubscribe from this group and all its topics, send an email to mosek+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/mosek/0ca111f6-385f-47c9-a23d-ec495b842394n%40googlegroups.com.

[Michal Adamaszek]

Write the model to an OPF file https://docs.mosek.com/9.2/cxxfusion/debugging-tutorials.html for small data and slowly check if it matches your expectations. 

Your conversion is at first sight correct, but a bit roundabout. You could just as well have done something like

t->index(0), 0.5, Expr::sub(x_0, start_)   belongs to   rotated quadratic cone

ie. why take a square root and then square it, when you can do it in one step.

If your investigation with the opf file does not help or if it raises questions then please let me know and time permitting I will try to look more carefully.

[G Prathap]

Hi Michal,

Thank you so much, I will have a look and let you know 

[Michal Adamaszek]

The main problems are:

1. You declare lambda values as int, so one gets rounded to 0.

2. You forget the connection between the 0-th and 1-st point, so there can be a jump directly from start to end. So the t variable should be longer by 1.

3. The order of some data inputs in the m file and cc file are different, so it is hard to compare exactly.

Of course 2. is the main issue.

Actually I would recommend writing the data to filename.ptf with the ptf extension. Then things are much more clearly visible.

[Michal Adamaszek]

And for the record here is a bit more Fusion-ish code for your problem, assuming that all norms are 2-norms (although in the m file you mix 1- and 2-norms; for 1-norm you would need a bit more complication).

#include <iostream>

#include "fusion.h"

using namespace mosek::fusion;

using namespace monty;

#define nint1(a)  new_array_ptr<int>({(a)})

#define nint(a,b) new_array_ptr<int>({(a),(b)})

int main(int argc, char ** argv)

{

  int number_of_steps = 15;

  double lambda1_val = 1000;

  double lambda2_val = 0.1;

  double lambda3_val = 1000;

  int dim_space = 3;

  auto A = new_array_ptr<double, 2>({{-0.470233,    0.882543,         0},

                                      {0.470233,   -0.882543,         0},

                                      {-0.807883,   -0.430453,    0.402535},

                                      {0.807883,    0.430453,   -0.402535},

                                      {-0.355254,   -0.189285,   -0.915405},

                                      {0.355254,    0.189285,    0.915405}});

 

  auto b = new_array_ptr<double, 1>({3.21407,

                                    0.785929,

                                    4.81961,

                                  -1.40824,

                                    0.878975,

                                    1.12103});

  auto end_ = new_array_ptr<double, 1>({-1, -1, 1});

  auto start_ = new_array_ptr<double, 1>({-4, 0, 2});

  Model::t M = new Model();

  auto x = M->variable("x", new_array_ptr<int,1>({dim_space, number_of_steps}), Domain::unbounded()) ;

  auto t      = M->variable("t", number_of_steps-1, Domain::unbounded());

  auto tstart = M->variable("ts",Domain::unbounded());

  auto tend   = M->variable("te",Domain::unbounded());

  M->constraint(Expr::vstack(tstart, 0.5, Expr::sub(x->slice(nint(0,0), nint(3,1))->reshape(nint1(3)),

                                                    start_)), 

                Domain::inRotatedQCone());

  M->constraint(Expr::hstack(t, 

                             Expr::constTerm(number_of_steps-1, 0.5), 

                             Expr::transpose(Expr::sub(x->slice(nint(0,0), nint(3,number_of_steps-1)),

                                                       x->slice(nint(0,1), nint(3,number_of_steps))))),

                Domain::inRotatedQCone());

  M->constraint(Expr::vstack(tend, 0.5, Expr::sub(x->slice(nint(0,number_of_steps-1), nint(3,number_of_steps))->reshape(nint1(3)),

                                                  end_)), 

                Domain::inRotatedQCone());

  for (int i = 0; i < number_of_steps; i++)

    M->constraint(Expr::mul(A, x->slice(nint(0,i), nint(3,i+1))->reshape(nint1(3))), Domain::lessThan(b));

  M->setLogHandler([](const std::string& msg){std::cout<< msg << std::flush;});

  auto lambda1 = M->parameter("lambda_1");

  auto lambda2 = M->parameter("lambda_2");

  auto lambda3 = M->parameter("lambda_3");

  lambda1->setValue(lambda1_val);

  lambda2->setValue(lambda2_val);

  lambda3->setValue(lambda3_val);

  M->objective(ObjectiveSense::Minimize, 

               Expr::add( Expr::sum(Expr::mul(lambda2, t)),

                          Expr::add(Expr::mul(lambda1, tstart), Expr::mul(lambda3, tend))));

  M->writeTask("a.ptf");

  M->solve();

  auto sol = *(x->level());

  std::cout<< "solution "<< sol << std::endl;

}

[G Prathap]

Thank you so much for your help. i really appreciate your effort 🙂👌 

To view this discussion on the web visit https://groups.google.com/d/msgid/mosek/cdc6675c-e106-4289-9dd5-9c396efde7f6n%40googlegroups.com.

[G Prathap]

First of all, thank you so much for the help so far, I want to use 1-norm 1 instead of 2-norm,

Matlab code using 1-norm 

clear all
close all
clc

number_of_steps = 5;
lambda1 = 1.0;
lambda2 = 1.0;
lambda3 = 1.0;

dim_ = 3;
Ab = [-0.470233  0.882543         0   3.21407
 0.470233 -0.882543        -0  0.785929
-0.807883 -0.430453  0.402535   4.81961
 0.807883  0.430453 -0.402535  -1.40824
-0.355254 -0.189285 -0.915405  0.878975
 0.355254  0.189285  0.915405   1.12103];

A = Ab(:,1:dim_);
b = Ab(:,dim_+1);


start = [-4 , 0, 2];
goal = [-1, -1, 1];

cvx_solver mosek
cvx_begin
    variable x(3, number_of_steps)

    cost = 0;
    for i = 1:(number_of_steps-1)

        cost = cost + norm(x(:, i) - x(:, i+1), 1)*lambda2;
    end
    cost = cost + norm(x(:, 1) - start', 1)*lambda1;
    cost = cost + norm(x(:, number_of_steps) - goal', 1)*lambda3;

    minimize(cost)
    subject to
        for i = 1:number_of_steps
            A*x(:, i) <= b;
        end
cvx_end

scatter3(x(1,:), x(2,:), x(3,:), 50, 'k', 'filled'); hold on

C++ Mosek using 1-norm 

    #include <iostream>
    #include "fusion.h"
    #include "monty.h"

    using namespace mosek::fusion;
    using namespace monty;

    #define nint1(a)  new_array_ptr<int>({(a)})

    #define nint(a, b) new_array_ptr<int>({(a),(b)})

    int main(int argc, char ** argv)
    {

      int num_steps = 5;

      Model::t M = new Model();

      int dim_space = 3;
      double lambda1_val = 1.0;
      double lambda2_val = 1.0;
      double lambda3_val = 1.0;

      auto A = new_array_ptr<double, 2>({{-0.470233,    0.882543,         0},
                                          {0.470233,   -0.882543,         0},
                                          {-0.807883,   -0.430453,    0.402535},
                                          {0.807883,    0.430453,   -0.402535},
                                          {-0.355254,   -0.189285,   -0.915405},
                                          {0.355254,    0.189285,    0.915405}});
     
      auto b = new_array_ptr<double, 1>({3.21407,
                                        0.785929,
                                        4.81961,
                                      -1.40824,
                                        0.878975,
                                        1.12103});

      auto end_ = new_array_ptr<double, 1>({-1, -1, 1});
      auto start_ = new_array_ptr<double, 1>({-4, 0, 2});

      auto x = M->variable("x", new_array_ptr<int,1>({dim_space, num_steps}), Domain::unbounded());
      auto z = M->variable("z", new_array_ptr<int,1>({dim_space, num_steps-1}), Domain::unbounded());
      auto z_norm = M->variable("z_norm", num_steps-1,Domain::unbounded());

      auto tstart = M->variable("ts", dim_space, Domain::unbounded());
      auto tstart_t = M->variable("ts_norm", Domain::unbounded());
      auto tend   = M->variable("te", dim_space, Domain::unbounded());
      auto tend_t   = M->variable("te_norm",Domain::unbounded());

      // Absolute value
      // t >= |x|, where t, x have the same shape
      // 1-norm
      // t >= sum( |x_i| ), x is a vector variable
     
      M->constraint(Expr::add(tstart,
                    Expr::sub(x->slice(nint(0,0), nint(dim_space, 1))->reshape(nint1(dim_space)), start_)), Domain::greaterThan(0.0));
      M->constraint(Expr::sub(tstart,
                    Expr::sub(x->slice(nint(0,0), nint(dim_space, 1))->reshape(nint1(dim_space)), start_)), Domain::greaterThan(0.0));
      M->constraint(Expr::sub(tstart_t, Expr::sum(tstart)), Domain::greaterThan(0.0));
     
      M->constraint(Expr::add(tend,
                    Expr::sub(x->slice(nint(0,num_steps-1), nint(dim_space,num_steps))->reshape(nint1(dim_space)), end_)), Domain::greaterThan(0.0));
      M->constraint(Expr::sub(tend,
                    Expr::sub(x->slice(nint(0,num_steps-1), nint(dim_space,num_steps))->reshape(nint1(dim_space)), end_)), Domain::greaterThan(0.0));
      M->constraint(Expr::sub(tend_t, Expr::sum(tend)), Domain::greaterThan(0.0));


      M->constraint(Expr::add(z, Expr::sub(x->slice(nint(0,0), nint(dim_space,num_steps-1)),x->slice(nint(0,1), nint(dim_space, num_steps))))
                                                            , Domain::greaterThan(0.0));
      M->constraint(Expr::sub(z, Expr::sub(x->slice(nint(0,0), nint(dim_space,num_steps-1)), x->slice(nint(0,1), nint(dim_space, num_steps))))
                                                            , Domain::greaterThan(0.0));
      M->constraint(Expr::sub(z_norm, Expr::sum(z, 0)), Domain::greaterThan(0.0));  

      for (int i = 0; i < num_steps; i++){
        M->constraint(Expr::mul(A, x->slice(nint(0,i), nint(3,i+1))->reshape(nint1(dim_space))), Domain::lessThan(b));
      }
     
      M->objective(ObjectiveSense::Minimize, Expr::add(Expr::sum(Expr::mul(lambda2_val, z_norm)), Expr::add(Expr::mul(lambda1_val, tstart_t), Expr::mul(lambda3_val, tend_t))));
      M->writeTask("problem.ptf");
      M->setLogHandler([](const std::string & msg) { std::cout << msg << std::flush; } );
      M->solve();
     
      std::cout<< *(x->level()) << std::endl;

    }

Output 

where black dots depict the result of Matlab implementation, C++ version all the time, i.e., for different step sizes, get only three points, i.e., start, end, one point with the constraint space (Ax<=b). 

If someone could explain to me where I did wrong in the 1-norm implementation. Here I attached the ptf file; it seems to be fine when I go through it. 

Thanks,

Geesara

[Michal Adamaszek]

There is nothing wrong. With the 1-norm the solution is highly non-unique, but the one from cvx and c++ are just as optimal (note that you don't t have three points, but some points repeat).

[G Prathap]

Sure, thank you 

To view this discussion on the web visit https://groups.google.com/d/msgid/mosek/da0fe86b-fd81-4aa6-9671-0f347e5494c8n%40googlegroups.com.