1 // 2 // Test.cpp 3 // 4 // 5 // Created by Hassan on 3 Jan 2016. 6 // 7 #define DOCTEST_CONFIG_IMPLEMENT_WITH_MAIN 8 #include <stdio.h> 9 #include <iostream> 10 #include <string> 11 #include <functional> 12 #include <stdio.h> 13 #include <cstring> 14 #include <fstream> 15 #include <stdio.h> 16 #include <stdlib.h> 17 #include <gravity/solver.h> 18 #include <gravity/doctest.h> 19 #include <PowerNet.h> 20 21 22 using namespace std; 23 using namespace gravity; 24 25 #ifdef USE_MP 26 TEST_CASE("testing readNL() function on ex4.nl") { 27 Model<> M; 28 string NL_file = string(prj_dir)+"/data_sets/NL/ex4.nl"; 29 int status = M.readNL(NL_file); 30 M.print(); 31 CHECK(status==0); 32 CHECK(M.get_nb_vars()==36); 33 CHECK(M.get_nb_cons()==30); 34 CHECK(M.is_convex()); 35 M.restructure(); 36 DebugOn("Done Restructuring!"); 37 M.write(); 38 } 39 40 TEST_CASE("testing readNL() function on waterlund32.nl") { 41 Model<> M; 42 string NL_file = string(prj_dir)+"/data_sets/NL/waterund32.nl"; 43 int status = M.readNL(NL_file); 44 M.print(); 45 CHECK(status==0); 46 CHECK(M.get_nb_vars()==660); 47 CHECK(M.get_nb_cons()==200); 48 CHECK(!M.is_convex()); 49 M.restructure(); 50 DebugOn("Done Restructuring!"); 51 } 52 #endif 53 54 TEST_CASE("testing param, var anf func copy operators") { 55 indices ids("ids"); 56 ids.add("key1", "key2"); 57 param<int> ip("ip"); 58 ip.in(ids); 59 ip.print(); 60 ip.set_val("key1",2); 61 ip.set_val("key2",3); 62 ip.print(); 63 CHECK(ip._range->first==2); 64 CHECK(ip._range->second==3); 65 param<int> ip2 = ip("key1");/* ip2 is a masked copy of ip, they will be sharing the same _val vector */ 66 ip2.print(); 67 CHECK(ip2.is_param()); 68 CHECK(ip2.get_intype()==integer_); 69 CHECK(ip2.is_positive()); 70 param<int> ip3 = ip.deep_copy();/* ip3 will not be sharing the same _val vector with ip or ip2 */ 71 ip.set_val("key1",-1); 72 CHECK(ip._range->first==-1);/* The range of ip should be updated */ 73 CHECK(ip2._range->first==-1);/* So is the range of ip2 */ 74 CHECK(ip3.eval("key1")==2);/* The _val of ip3 was not modified */ 75 var<> v1("v1", -2,3), v2("v2", -10,40); 76 v1.in(ids);v2.in(ids); 77 func<> f = v1 + v1.get_lb()*v2; 78 f.print(); 79 auto fcpy = f; 80 func<> deep_cpy; 81 deep_cpy.deep_copy(f); 82 fcpy.print(); 83 v1.set_lb("key1", -3); 84 f.uneval(); 85 f.eval_all(); 86 fcpy.uneval(); 87 fcpy.eval_all(); 88 f.print(); 89 fcpy.print(); 90 deep_cpy.print(); 91 CHECK(fcpy.to_str(0,5)=="v1[key1] - 3v2[key1]"); 92 CHECK(deep_cpy.to_str(0,5)=="v1[key1] - 2v2[key1]"); 93 } 94 95 96 //TEST_CASE("testing projection3") { 97 // string fname = string(prj_dir)+"/data_sets/Power/pglib_opf_case3_lmbd.m"; 98 //// string fname = string(prj_dir)+"/data_sets/Power/nesta_case5_pjm.m"; 99 // int output = 0; 100 // double tol = 1e-6; 101 // PowerNet grid; 102 // grid.readgrid(fname); 103 // auto ACOPF = build_ACOPF(grid,ACRECT); 104 // ACOPF->print_symbolic(); 105 // ACOPF->print(); 106 // 107 // ACOPF->project(); 108 // ACOPF->print(); 109 //// CHECK(Mtest.get_nb_cons() == 4); 110 //} 111 112 //m = Model(solver=solver) 113 // 114 //@variable(m, x[1:8]) 115 // 116 //setlowerbound(x[1], 100) 117 //setlowerbound(x[2], 1000) 118 //setlowerbound(x[3], 1000) 119 //setlowerbound(x[4], 10) 120 //setlowerbound(x[5], 10) 121 //setlowerbound(x[6], 10) 122 //setlowerbound(x[7], 10) 123 //setlowerbound(x[8], 10) 124 // 125 //setupperbound(x[1], 10000) 126 //setupperbound(x[2], 10000) 127 //setupperbound(x[3], 10000) 128 //setupperbound(x[4], 1000) 129 //setupperbound(x[5], 1000) 130 //setupperbound(x[6], 1000) 131 //setupperbound(x[7], 1000) 132 //setupperbound(x[8], 1000) 133 // 134 //@constraint(m, 0.0025*(x[4] + x[6]) <= 1) 135 //@constraint(m, 0.0025*(x[5] - x[4] + x[7]) <= 1) 136 //@constraint(m, 0.01(x[8]-x[5]) <= 1) 137 //@NLconstraint(m, 100*x[1] - x[1]*x[6] + 833.33252*x[4] <= 83333.333) 138 //@NLconstraint(m, x[2]*x[4] - x[2]*x[7] - 1250*x[4] + 1250*x[5] <= 0) 139 //@NLconstraint(m, x[3]*x[5] - x[3]*x[8] - 2500*x[5] + 1250000 <= 0) 140 // 141 //@objective(m, Min, x[1]+x[2]+x[3]) 142 143 144 TEST_CASE("Variable Scaling") { 145 Model<> M("Test"); 146 param<> lb("x_lb"); 147 lb = {100,1000,1000,10,10,10,10,10}; 148 param<> ub("x_lb"); 149 ub = {10000,10000,10000,1000,1000,1000,1000,1000}; 150 var<> x("x",lb,ub); 151 M.add(x.in(range(1,8))); 152 153 Constraint<> C1("C1"); 154 C1 = 0.0025*(x[4] + x[6]); 155 M.add(C1 <= 1); 156 157 Constraint<> C2("C2"); 158 C2 = 0.0025*(x[5] - x[4] + x[7]); 159 M.add(C2 <= 1); 160 161 Constraint<> C3("C3"); 162 C3 = 0.01*(x[8]-x[5]); 163 M.add(C3 <= 1); 164 165 Constraint<> C4("C4"); 166 C4 = 100*x[1] - x[1]*x[6] + 833.33252*x[4]; 167 M.add(C4 <= 83333.333); 168 Constraint<> C5("C5"); 169 C5 = x[2]*x[4] - x[2]*x[7] - 1250*x[4] + 1250*x[5]; 170 M.add(C5 <= 0); 171 Constraint<> C6("C6"); 172 C6 = x[3]*x[5] - x[3]*x[8] - 2500*x[5] + 1250000; 173 M.add(C6 <= 0); 174 175 M.min(x[1]+x[2]+x[3]); 176 M.print(); 177 178 solver<> s(M,ipopt); 179 s.run(5,1e-6); 180 auto obj_val = M.get_obj_val(); 181 auto M_scale = M; 182 M_scale.reset(); 183 M_scale.scale_vars(10); 184 M_scale.print(); 185 double coef_scale = 100; 186 M_scale.scale_coefs(coef_scale); 187 M_scale.print(); 188 M_scale.initialize_midpoint(); 189 solver<> s2(M_scale,ipopt); 190 s2.run(5,1e-6); 191 CHECK(std::abs(M_scale.get_obj_val()-obj_val) < 1e-3); 192 M_scale.write_solution(); 193 } 194 195 TEST_CASE("Model.relax()") { 196 Model<> M("Test"); 197 param<> lb("x_lb"); 198 lb = {100,1000,1000,10,10,10,10,10}; 199 param<> ub("x_lb"); 200 ub = {10000,10000,10000,1000,1000,1000,1000,1000}; 201 var<> x("x",lb,ub); 202 auto x_ids = indices("x_ids"); 203 x_ids = range(1,8); 204 M.add(x.in(x_ids)); 205 206 indices x_mat("x_mat"), x_mat_RLT1("x_mat_RLT1"), x_mat_RLT2("x_mat_RLT2"), x_mat_RLT3("x_mat_RLT3"); 207 208 param<> A("A"); 209 A.in(x_ids); 210 A._indices->add_in_row(0, "1"); 211 A._indices->add_in_row(0, "2"); 212 x_mat.add_in_row(0, "4"); 213 x_mat.add_in_row(0, "6"); 214 x_mat_RLT1.add_in_row(0, "1"); 215 x_mat_RLT1.add_in_row(0, "1"); 216 x_mat_RLT2.add_in_row(0, "2"); 217 x_mat_RLT2.add_in_row(0, "2"); 218 x_mat_RLT3.add_in_row(0, "3"); 219 x_mat_RLT3.add_in_row(0, "3"); 220 A.set_val("1", 0.0025);A.set_val("2", 0.0025); 221 A._indices->add_in_row(1, "3"); 222 A._indices->add_in_row(1, "4"); 223 A._indices->add_in_row(1, "5"); 224 x_mat.add_in_row(1, "5"); 225 x_mat.add_in_row(1, "4"); 226 x_mat.add_in_row(1, "7"); 227 x_mat_RLT1.add_in_row(1, "1"); 228 x_mat_RLT1.add_in_row(1, "1"); 229 x_mat_RLT1.add_in_row(1, "1"); 230 x_mat_RLT2.add_in_row(1, "2"); 231 x_mat_RLT2.add_in_row(1, "2"); 232 x_mat_RLT2.add_in_row(1, "2"); 233 x_mat_RLT3.add_in_row(1, "3"); 234 x_mat_RLT3.add_in_row(1, "3"); 235 x_mat_RLT3.add_in_row(1, "3"); 236 A.set_val("3", 0.0025);A.set_val("4", -0.0025);A.set_val("5", 0.0025); 237 A._indices->add_in_row(2, "6"); 238 A._indices->add_in_row(2, "7"); 239 x_mat.add_in_row(2, "8"); 240 x_mat.add_in_row(2, "5"); 241 x_mat_RLT1.add_in_row(2, "1"); 242 x_mat_RLT1.add_in_row(2, "1"); 243 x_mat_RLT2.add_in_row(2, "2"); 244 x_mat_RLT2.add_in_row(2, "2"); 245 x_mat_RLT3.add_in_row(2, "3"); 246 x_mat_RLT3.add_in_row(2, "3"); 247 A.set_val("6", 0.01);A.set_val("7", -0.01); 248 249 Constraint<> LinCons("LinCons"); 250 LinCons = A*x.in(x_mat); 251 M.add(LinCons.in(range(1,3)) <= 1); 252 253 Constraint<> RLT1("RLT1"); 254 RLT1 = A*x.in(x_mat_RLT1)*x.in(x_mat) - x.in(x.repeat_id(3,0)); 255 M.add(RLT1.in(range(1,3)) <= 0); 256 257 Constraint<> RLT2("RLT2"); 258 RLT2 = A*x.in(x_mat_RLT2)*x.in(x_mat) - x.in(x.repeat_id(3,1)); 259 M.add(RLT2.in(range(1,3)) <= 0); 260 261 Constraint<> RLT3("RLT3"); 262 RLT3 = A*x.in(x_mat_RLT3)*x.in(x_mat) - x.in(x.repeat_id(3,2)); 263 M.add(RLT3.in(range(1,3)) <= 0); 264 265 /* 266 Constraint<> C1("C1"); 267 C1 = 0.0025*(x[4] + x[6]); 268 M.add(C1 <= 1); 269 270 Constraint<> C2("C2"); 271 C2 = 0.0025*(x[5] - x[4] + x[7]); 272 M.add(C2 <= 1); 273 274 Constraint<> C3("C3"); 275 C3 = 0.01*(x[8]-x[5]); 276 M.add(C3 <= 1); 277 */ 278 279 280 Constraint<> C4("C4"); 281 C4 = 100*x[1] - x[1]*x[6] + 833.33252*x[4]; 282 M.add(C4 <= 83333.333); 283 Constraint<> C5("C5"); 284 C5 = x[2]*x[4] - x[2]*x[7] - 1250*x[4] + 1250*x[5]; 285 M.add(C5 <= 0); 286 Constraint<> C6("C6"); 287 C6 = x[3]*x[5] - x[3]*x[8] - 2500*x[5] + 1250000; 288 M.add(C6 <= 0); 289 290 291 M.min(x[1]+x[2]+x[3]); 292 293 // M.scale_vars(100); 294 // double coef_scale = 100; 295 // M.scale_coefs(coef_scale); 296 M.print(); 297 298 auto determinant_level = 1; 299 bool add_Kim_Kojima = false, add_SDP_3d = false; 300 auto LB = M.relax(determinant_level,add_Kim_Kojima, add_SDP_3d); 301 302 // LB->print(); 303 // LB->scale_vars(100); 304 // LB->print(); 305 // double coef_scale = 100; 306 // LB->scale_coefs(coef_scale); 307 LB->print(); 308 // M.print(); 309 double max_time = 54000,ub_solver_tol=1e-6, lb_solver_tol=1e-6, range_tol=1e-4, opt_rel_tol=1e-2, opt_abs_tol=1e6; 310 unsigned max_iter=30, nb_threads = thread::hardware_concurrency(); 311 SolverType ub_solver_type = ipopt, lb_solver_type = ipopt; 312 M.run_obbt(LB, max_time, max_iter, opt_rel_tol, opt_abs_tol, nb_threads=1, ub_solver_type, lb_solver_type, ub_solver_tol, lb_solver_tol, range_tol); 313 LB->print_constraints_stats(1e-6); 314 // LB->print_nonzero_constraints(1e-6); 315 LB->print(); 316 auto final_gap = 100*(M.get_obj_val() - LB->get_obj_val())/std::abs(M.get_obj_val()); 317 CHECK(final_gap<1); 318 } 319 320 TEST_CASE("hard nlp") { 321 Model<> M("Test"); 322 param<> lb("x_lb"); 323 lb = {100,1000,1000,10,10,10,10,10}; 324 param<> ub("x_lb"); 325 ub = {10000,10000,10000,1000,1000,1000,1000,1000}; 326 var<> x("x",lb,ub); 327 M.add(x.in(range(1,8))); 328 329 330 Constraint<> C1("C1"); 331 C1 = 0.0025*(x[4] + x[6]); 332 M.add(C1 <= 1); 333 334 335 Constraint<> RLT1_1("RLT1_1"); 336 RLT1_1 = x[1]*(0.0025*(x[4] + x[6]) - 1); 337 M.add(RLT1_1 <= 0); 338 339 Constraint<> RLT1_2("RLT1_2"); 340 RLT1_2 = x[2]*(0.0025*(x[4] + x[6]) - 1); 341 M.add(RLT1_2 <= 0); 342 343 Constraint<> RLT1_3("RLT1_3"); 344 RLT1_3 = x[3]*(0.0025*(x[4] + x[6]) - 1); 345 M.add(RLT1_3 <= 0); 346 347 348 Constraint<> C2("C2"); 349 C2 = 0.0025*(x[5] - x[4] + x[7]); 350 M.add(C2 <= 1); 351 352 353 Constraint<> RLT2_1("RLT2_1"); 354 RLT2_1 = x[1]*(0.0025*(x[5] - x[4] + x[7])-1); 355 M.add(RLT2_1 <= 0); 356 357 Constraint<> RLT2_2("RLT2_2"); 358 RLT2_2 = x[2]*(0.0025*(x[5] - x[4] + x[7])-1); 359 M.add(RLT2_2 <= 0); 360 361 Constraint<> RLT2_3("RLT2_3"); 362 RLT2_3 = x[3]*(0.0025*(x[5] - x[4] + x[7])-1); 363 M.add(RLT2_3 <= 0); 364 365 366 Constraint<> C3("C3"); 367 C3 = 0.01*(x[8]-x[5]); 368 M.add(C3 <= 1); 369 370 371 Constraint<> RLT3_1("RLT3_1"); 372 RLT3_1 = x[1]*(0.01*(x[8]-x[5])-1); 373 M.add(RLT3_1 <= 0); 374 375 Constraint<> RLT3_2("RLT3_2"); 376 RLT3_2 = x[2]*(0.01*(x[8]-x[5])-1); 377 M.add(RLT3_2 <= 0); 378 379 Constraint<> RLT3_3("RLT3_3"); 380 RLT3_3 = x[3]*(0.01*(x[8]-x[5])-1); 381 M.add(RLT3_3 <= 0); 382 383 Constraint<> C4("C4"); 384 C4 = 100*x[1] - x[1]*x[6] + 833.33252*x[4]; 385 M.add(C4 <= 83333.333); 386 Constraint<> C5("C5"); 387 C5 = x[2]*x[4] - x[2]*x[7] - 1250*x[4] + 1250*x[5]; 388 M.add(C5 <= 0); 389 Constraint<> C6("C6"); 390 C6 = x[3]*x[5] - x[3]*x[8] - 2500*x[5] + 1250000; 391 M.add(C6 <= 0); 392 393 394 M.min(x[1]+x[2]+x[3]); 395 396 // M.scale_vars(100); 397 // double coef_scale = 100; 398 // M.scale_coefs(coef_scale); 399 M.print(); 400 401 auto determinant_level = 1; 402 bool add_Kim_Kojima = false, add_SDP_3d = false; 403 auto LB = M.relax(determinant_level,add_Kim_Kojima, add_SDP_3d); 404 405 // LB->print(); 406 // LB->scale_vars(100); 407 // LB->print(); 408 // double coef_scale = 100; 409 // LB->scale_coefs(coef_scale); 410 LB->print(); 411 // M.print(); 412 double max_time = 54000,ub_solver_tol=1e-6, lb_solver_tol=1e-6, range_tol=1e-4, opt_rel_tol=1e-2, opt_abs_tol=1e6; 413 unsigned max_iter=1, nb_threads = thread::hardware_concurrency(); 414 SolverType ub_solver_type = ipopt, lb_solver_type = ipopt; 415 M.run_obbt(LB, max_time, max_iter, opt_rel_tol, opt_abs_tol, nb_threads=1, ub_solver_type, lb_solver_type, ub_solver_tol, lb_solver_tol, range_tol); 416 LB->print_constraints_stats(1e-6); 417 // LB->print_nonzero_constraints(1e-6); 418 LB->print(); 419 } 420 421 TEST_CASE("testing projection1") { 422 indices buses("buses"); 423 buses.insert("1", "2", "3", "4"); 424 indices node_pairs("bpairs"); 425 node_pairs.insert("1,2", "1,3", "3,4", "4,1"); 426 427 Model<> Mtest("Mtest"); 428 var<> R_Wij("R_Wij", -1, 1); 429 /* Imaginary part of Wij = ViVj */ 430 var<> Im_Wij("Im_Wij", -1, 1); 431 var<> Wii("Wii", 0.8, 1.21); 432 Mtest.add(R_Wij.in(node_pairs), Im_Wij.in(node_pairs), Wii.in(buses)); 433 Constraint<> SOC("SOC"); 434 SOC = 2*R_Wij + pow(Im_Wij, 2) - 4*Wii.from(node_pairs); 435 Mtest.add(SOC.in(node_pairs) == 0); 436 437 auto subset = node_pairs.exclude("4,1"); 438 Constraint<> PAD("PAD"); 439 PAD = 2*R_Wij.in(subset) - Im_Wij.in(subset); 440 Mtest.add(PAD.in(subset)<=2); 441 442 443 Constraint<> PAD2("PAD2"); 444 PAD2 = R_Wij("4,1") - 2*Im_Wij("4,1"); 445 Mtest.add(PAD2>=1); 446 447 Mtest.min(sum(R_Wij)); 448 Mtest.print(); 449 CHECK(Mtest.get_nb_cons() == 8); 450 Mtest.project(); 451 Mtest.print(); 452 CHECK(Mtest.get_nb_eq() == 0); 453 CHECK(Mtest.get_nb_ineq() == 12); 454 } 455 456 TEST_CASE("testing projection2") { 457 indices buses("buses"); 458 buses.insert("1", "2", "3", "4"); 459 indices node_pairs("bpairs"); 460 node_pairs.insert("1,2", "1,3", "3,4", "4,1"); 461 auto subset = node_pairs.exclude("4,1"); 462 463 Model<> Mtest("Mtest"); 464 var<> R_Wij("R_Wij"); 465 /* Imaginary part of Wij = ViVj */ 466 var<> Im_Wij("Im_Wij", -1, 1); 467 var<> Wii("Wii", 0.8, 1.21); 468 Mtest.add(R_Wij.in(node_pairs), Im_Wij.in(node_pairs), Wii.in(buses)); 469 470 Constraint<> SOC("SOC"); 471 SOC = 2*R_Wij.in(subset) + pow(Im_Wij.in(subset), 2) - 4*Wii.from(subset); 472 Mtest.add(SOC.in(subset) == 0); 473 474 475 Constraint<> PAD("PAD"); 476 PAD = 2*R_Wij.in(node_pairs) - Im_Wij.in(node_pairs); 477 Mtest.add(PAD.in(node_pairs)<=2); 478 479 480 Mtest.min(sum(R_Wij)); 481 Mtest.print(); 482 CHECK(Mtest.get_nb_cons() == 7); 483 Mtest.project(); 484 Mtest.print(); 485 CHECK(Mtest.get_nb_eq() == 0); 486 CHECK(Mtest.get_nb_ineq() == 4); 487 } 488 // 489 // 490 // 491 //TEST_CASE("testing projection3") { 492 // string fname = string(prj_dir)+"/data_sets/Power/pglib_opf_case3_lmbd.m"; 493 // int output = 0; 494 // double tol = 1e-6; 495 // PowerNet grid; 496 // grid.readgrid(fname); 497 // auto SOCOPF = grid.build_SCOPF(); 498 // SOCOPF->print_symbolic(); 499 // SOCOPF->print(); 500 // 501 // SOCOPF->project(); 502 // SOCOPF->print(); 503 //// CHECK(Mtest.get_nb_cons() == 4); 504 //} 505 506 507 TEST_CASE("testing constants") { 508 constant<> c0; 509 CHECK(c0.is_double()); 510 c0 = 3.5; 511 CHECK(c0==3.5); 512 CHECK(c0.is_positive()); 513 CHECK(!c0.is_negative()); 514 constant<Cpx> cx0; 515 cx0 = Cpx(-1,1); 516 CHECK(cx0.is_complex()); 517 CHECK(cx0==Cpx(-1,1)); 518 constant<Cpx> cx1; 519 cx1 = Cpx(-1,-2); 520 auto cx2 = cx0 + cx1; 521 cx2.print(); 522 CHECK(cx2.is_complex()); 523 CHECK(cx2.is_negative()); 524 auto mag0 = sqrmag(cx2); 525 CHECK(std::abs(mag0.eval()-5)<1e-8); 526 auto ang0 = angle(cx2); 527 ang0.println(); 528 CHECK(std::abs(ang0.eval()-(-2.677945045))<1e-8); 529 auto cx3 = conj(cx1); 530 CHECK(cx3.eval().real()==cx1.eval().real()); 531 CHECK(cx3.eval().imag()==-1*cx1.eval().imag()); 532 CHECK(real(cx3) == cx3.eval().real()); 533 CHECK(imag(cx3) == cx3.eval().imag()); 534 CHECK(cx3.get_dim() == 1); 535 try{ 536 cx3.get_dim(2); 537 } 538 catch(invalid_argument& arg){ 539 cout << "Error successfully caught: "<< endl; 540 cout << arg.what() << endl; 541 } 542 CHECK(cx3.is_number()); 543 } 544 545 TEST_CASE("testing parameters") { 546 param<> c0("c0"); 547 CHECK(c0.is_double()); 548 c0 = 3.5; 549 CHECK(c0.is_positive()); 550 CHECK(!c0.is_negative()); 551 param<Cpx> cx0("cx0"); 552 cx0 = Cpx(-1,1); 553 CHECK(cx0.is_complex()); 554 param<Cpx> cx1("cx1"); 555 cx1 = Cpx(-1,-2); 556 auto mag0 = sqrmag(cx1); 557 auto ang0 = ang(cx1); 558 mag0.print(); 559 ang0.print(); 560 auto cx2 = conj(cx1); 561 CHECK(cx2.eval().real()==cx1.eval().real()); 562 CHECK(cx2.get_dim() == 1); 563 param<Cpx> cx3("cx3"); 564 cx3.in(C(6)); 565 CHECK(cx3.get_dim() == 6); 566 } 567 568 TEST_CASE("testing param copy operator") { 569 param<int> ip("ip"); 570 ip.print(); 571 ip.add_val(2); 572 ip.add_val(3); 573 ip.print(); 574 param<int> ip2(ip); 575 ip2.print(); 576 CHECK(ip==ip2); 577 CHECK(ip2.is_param()); 578 CHECK(ip2.get_intype()==integer_); 579 CHECK(ip2.is_positive()); 580 CHECK(ip2._range->first==2); 581 CHECK(ip2._range->second==3); 582 } 583 584 TEST_CASE("testing param indexing, add_val() and set_val() functions") { 585 param<> ip("ip"); 586 ip.print(); 587 ip.add_val(2); 588 ip.add_val(-1.3); 589 CHECK(ip._range->first==-1.3); 590 CHECK(ip._range->second==2); 591 ip.print(); 592 ip.set_val(1, 1.5); 593 ip.print(); 594 CHECK(ip.eval(1)==1.5); 595 CHECK(ip._range->first==1.5); 596 CHECK(ip._range->second==2); 597 indices ids("index_set"); 598 ids.add({"id1", "id2", "key3"}); 599 param<> dp("dp"); 600 dp.in(ids); 601 dp.print(); 602 dp("id1") = 1.5; 603 dp("id2") = -231.5; 604 dp.print(); 605 CHECK(dp.eval("id1")==1.5); 606 CHECK(dp.eval("id2")==-231.5); 607 CHECK(dp.eval("key3")==0); 608 CHECK(dp._range->first==-231.5); 609 CHECK(dp._range->second==1.5); 610 REQUIRE_THROWS_AS(dp("unexisting_key").eval(), invalid_argument); 611 auto ndp = dp.in(ids.exclude("id2")); 612 ndp.print(); 613 CHECK(ndp.get_dim()==2); 614 } 615 616 TEST_CASE("testing matrix params") { 617 param<> mat("M"); 618 mat.set_size(3,4); 619 for (auto i = 0; i<3;i++) { 620 for (auto j = 0; j<4;j++) { 621 mat.set_val(i, j, 10*i+j); 622 } 623 } 624 mat.print(); 625 CHECK(mat.eval(1,2)==12); 626 CHECK(mat(1,2).eval()==12); 627 auto tr_mat = mat.tr(); 628 tr_mat.print(); 629 CHECK(tr_mat.eval(2,1)==12); 630 CHECK(tr_mat(2,1).eval()==12); 631 632 var<> v("v",0,1); 633 v.in(R(4)); 634 Constraint<> Mv("Mv"); 635 Mv = product(mat,v); 636 Mv.print(); 637 638 /* Complex matrices */ 639 param<Cpx> Cmat("Cmat"); 640 Cmat.set_size(3,3); 641 Cpx cval = Cpx(-1,2); 642 Cmat.set_val(0, 1, cval); 643 Cmat.print(); 644 CHECK(Cmat.eval(0,1)==cval); 645 CHECK(Cmat(0,1).eval()==cval); 646 } 647 648 TEST_CASE("testing param sign functions") { 649 param<int> ip("ip"); 650 ip.print(); 651 ip.add_val(2); 652 ip.add_val(3); 653 ip.print(); 654 CHECK(ip.is_non_negative()); 655 CHECK(ip.is_positive()); 656 CHECK(!ip.is_negative()); 657 CHECK(!ip.is_non_positive()); 658 param<> dp("dp"); 659 dp.add_val(0); 660 dp.add_val(-3); 661 dp.print(); 662 CHECK(!dp.is_positive()); 663 CHECK(!dp.is_non_negative()); 664 CHECK(!dp.is_negative()); 665 CHECK(dp.is_non_positive()); 666 } 667 668 TEST_CASE("testing variables indexing") { 669 indices ids("index_set"); 670 ids.add({"id1", "id2", "key3"}); 671 var<> iv("iv",-2, 5); 672 iv.in(ids); 673 iv.print(); 674 param<Cpx> lb("lb"), ub("ub"); 675 lb.in(ids); 676 lb = Cpx(0,-1); 677 ub.in(ids); 678 ub = Cpx(1,1); 679 var<Cpx> cv("cv", lb, ub); 680 cv.in(ids.exclude("id2")); 681 cv.print(); 682 CHECK(cv.get_dim()==2); 683 var<Cpx> cv1("cv", Cpx(0,-1),Cpx(1,1)); 684 cv1.in(ids.exclude("id2")); 685 cv1.print(); 686 CHECK(cv1.get_dim()==2); 687 CHECK(cv==cv1); 688 CHECK(cv.get_indices()!=ids); 689 } 690 691 692 693 TEST_CASE("testing vector dot product"){ 694 param<> a("a"); 695 a.set_size(3); 696 a.set_val(0, 1); 697 a.set_val(1, -1); 698 a.set_val(2, 2); 699 param<int> b("b"); 700 b=5; 701 CHECK(b.get_dim()==1); 702 CHECK(b.is_integer()); 703 CHECK(b.eval(0)==5); 704 b=0; 705 b=2; 706 b=5; 707 b=4; 708 CHECK(b.get_dim()==5); 709 CHECK(b.eval(1)==0); 710 b.set_val(1, 3); 711 CHECK(b.eval(1)==3); 712 b.print(); 713 var<> z("z",-1,1); 714 z.in(R(4)); 715 CHECK(z.get_dim()==4); 716 z.print(); 717 z.in(R(3)); 718 CHECK(z.get_dim()==3); 719 z.print(); 720 param<Cpx> cp("cp"); 721 cp.in(C(3)); 722 cp.set_val(Cpx(1,1)); 723 auto lin = cp+z; 724 CHECK(!lin.is_evaluated()); 725 lin.print_symbolic(); 726 lin.print(); 727 CHECK(lin._range->first==Cpx(0,1)); 728 CHECK(lin._range->second==Cpx(2,1)); 729 CHECK(lin.is_complex()); 730 CHECK(lin.is_linear()); 731 auto lin2 = cp*z; 732 auto lin3 = a.tr()*z - b; 733 CHECK(lin3.is_double()); 734 CHECK(lin3.get_dim()==5); 735 lin3.print(); 736 var<Cpx> y("y"); 737 y.in(C(5)); 738 auto lin4 = b*y*y; 739 lin4.print_symbolic(); 740 lin4.print(); 741 CHECK(lin4.is_complex()); 742 CHECK(lin4.is_convex()); 743 CHECK(lin4.get_dim()==5); 744 auto lin5 = lin4 - lin3; 745 lin5.print_symbolic(); 746 lin5.print(); 747 CHECK(lin5.is_complex()); 748 CHECK(lin5.get_dim()==5); 749 CHECK(lin5.get_nb_vars()==4); 750 CHECK(lin5.is_convex()); 751 CHECK(lin5.to_str()=="(b)y² - ([a]ᵀ)[z] + b"); 752 auto lin6 = (2*a-exp(a)+1).tr()*z; 753 lin6.print_symbolic(); 754 CHECK(lin6.to_str()=="([2a + 1 + -exp(a)]ᵀ)[z]"); 755 CHECK(lin6.is_linear()); 756 lin6.print(); 757 } 758 759 TEST_CASE("testing complex numbers") { 760 indices ids("index_set"); 761 ids.add({"id1", "id2", "key3", "key4"}); 762 var<> iv("x",-2, 5); 763 iv.in(ids); 764 var<Cpx> cv("y", Cpx(0,-1),Cpx(1,1)); 765 constant<Cpx> cx(Cpx(-1,-2)); 766 auto cx_conj = conj(cx); 767 CHECK(real(cx_conj)==real(cx)); 768 CHECK(imag(cx_conj).eval()==-1*imag(cx).eval()); 769 param<Cpx> px("px"); 770 px = Cpx(-1,-2); 771 px.print(); 772 auto px_conj = conj(px); 773 CHECK(real(px_conj)._is_real); 774 CHECK(imag(px_conj)._is_imag); 775 } 776 777 TEST_CASE("testing complex functions") { 778 var<> z("z",-1,1); 779 z.in(R(3)); 780 param<Cpx> cpx("cpx"); 781 cpx = Cpx(1,1); 782 cpx = Cpx(2,2); 783 cpx = Cpx(3,1); 784 auto cpx_f = 2*exp(cpx)*z; 785 cpx_f.print_symbolic(); 786 CHECK(cpx_f.to_str()=="((2,0)exp(cpx))z"); 787 cpx_f.print(); 788 } 789 790 TEST_CASE("testing ReLU") { 791 var<> A("A",-1,1), B("B",-1,1); 792 A.in(R(4)); 793 B.in(R(1)); 794 param<> X("X"); 795 X.in(R(4)); 796 X(0) = 0.2; 797 X(1) = 0.3; 798 X(2) = 0.4; 799 X(3) = -0.2; 800 auto f = ReLU(A.tr()*X + B); 801 f.print_symbolic(); 802 f.print(); 803 auto dfdA = f.get_derivative(B); 804 dfdA.print_symbolic(); 805 CHECK(dfdA.to_str()=="UnitStep(B + (Xᵀ)[A])"); 806 dfdA.print(); 807 } 808 809 TEST_CASE("2d Polynomial") { 810 var<> x1("x1",-1,1), x2("x2",-1,1); 811 Model<> M("test"); 812 M.add(x1.in(R(1)),x2.in(R(1))); 813 M.initialize_uniform(); 814 M.min(pow(x1,2) - 2*pow(x1,4) + x1*x2 - 4*pow(x2,2) + 4*pow(x2,8)); 815 solver<> s(M,ipopt); 816 s.set_option("max_iter", 2); 817 s.run(5,1e-6); 818 M.print(); 819 M.print_solution(); 820 } 821 822 TEST_CASE("testing range propagation") { 823 indices ids("index_set"); 824 ids.add({"id1", "id2", "key3", "key4"}); 825 var<> x("x",-2, 5); 826 x.in(ids); 827 auto f = 2*x + 2; 828 f.print_symbolic(); 829 CHECK(f.is_linear()); 830 CHECK(f.is_convex()); 831 CHECK(f._range->first==-2); 832 CHECK(f._range->second==12); 833 f+= pow(x,2); 834 CHECK(f._range->first==-2); 835 CHECK(f._range->second==37); 836 f+=2; 837 CHECK(f._range->first==0); 838 CHECK(f._range->second==39); 839 CHECK(f.to_str()=="x² + 2x + 4"); 840 f.print_symbolic(); 841 f.print(); 842 CHECK(f.is_quadratic()); 843 CHECK(f.is_convex()); 844 auto dfx = f.get_derivative(x); 845 dfx.print(); 846 CHECK(dfx.is_linear()); 847 CHECK(dfx.get_dim()==4); 848 param<> a("a"); 849 a.in(R(4)); 850 a = 2; 851 a.print(); 852 var<> x1("x1",-2, 2); 853 x1.in(ids); 854 var<> x2("x2",0, 3); 855 x2.in(ids); 856 auto f2 = a*pow(x1,4) * pow(x2,2); 857 f2.print_symbolic(); 858 f2.print(); 859 CHECK(f2.is_polynomial()); 860 CHECK(f2._range->first==0); 861 CHECK(f2._range->second==2*pow(2,4)*pow(3,2)); 862 var<> x3("x3",-2, 2); 863 x3.in(ids); 864 f2 *= pow(x3,3); 865 f2.print_symbolic(); 866 f2.print(); 867 CHECK(f2._range->first==-2*pow(2,4)*pow(3,2)*pow(2,3)); 868 CHECK(f2._range->second==2*pow(2,4)*pow(3,2)*pow(2,3)); 869 auto dfdx1 = f2.get_derivative(x1); 870 dfdx1.print_symbolic(); 871 dfdx1.print(); 872 CHECK(dfdx1._range->first==-8*pow(2,3)*pow(3,2)*pow(2,3)); 873 CHECK(dfdx1._range->second==8*pow(2,3)*pow(3,2)*pow(2,3)); 874 } 875 876 TEST_CASE("testing polynomial functions") { 877 param<int> a("a"); 878 a.set_size(3); 879 a.set_val(0,1); 880 a.set_val(1,-1); 881 a.set_val(2,2); 882 param<int> b("b"); 883 b=5; 884 b=0; 885 b=2; 886 var<> x("x",-1,1); 887 x.in(R(3)); 888 var<> y("y",-1,1); 889 y.in(R(3)); 890 var<> z("z",-1,1); 891 z.in(R(3)); 892 auto poly = pow(x,2)*pow(y,3)*pow(z,4) + b*pow(y,2)*pow(z,3); 893 CHECK(poly.to_str()=="x²y³z⁴ + (b)y²z³"); 894 poly.print_symbolic(); 895 poly.print(); 896 poly += a*x; 897 poly.print_symbolic(); 898 poly.print(); 899 auto dfdx = poly.get_derivative(x); 900 CHECK(dfdx.to_str()=="2xy³z⁴ + a"); 901 dfdx.print_symbolic(); 902 dfdx.print(); 903 auto dfd2x = dfdx.get_derivative(x); 904 CHECK(dfd2x.to_str()=="2y³z⁴"); 905 dfd2x.print_symbolic(); 906 dfd2x.print(); 907 } 908 909 910 TEST_CASE("testing complex matrix product") { 911 var<Cpx> X("X", Cpx(0,-1),Cpx(1,1)); 912 X.in(C(5,5)); 913 param<Cpx> A("A"); 914 A.set_size(3, 5); 915 A.set_val(0,0,Cpx(-1,1)); 916 A.set_val(1,0,Cpx(2,1)); 917 A.set_val(2,1,Cpx(-1,-1)); 918 A.set_val(1,1,Cpx(1,1)); 919 auto f = A*X; 920 f.print_symbolic(); 921 f.print(); 922 CHECK(f.is_linear()); 923 CHECK(f.is_convex()); 924 CHECK(f.get_dim()==15); 925 var<Cpx> Y("Y", Cpx(0,-1),Cpx(1,1)); 926 Y.in(C(5,5)); 927 f+= A*Y; 928 f.print_symbolic(); 929 f.print(); 930 CHECK(f.is_complex()); 931 CHECK(f.get_nb_vars()==50); 932 auto f2 = X*X; 933 f2.print_symbolic(); 934 f2.print(); 935 CHECK(f2.is_convex()); 936 f2 -= X*Y; 937 f2.print_symbolic(); 938 f2.print(); 939 CHECK(!f2.is_convex()); 940 CHECK(f2.get_dim()==25); 941 auto f3 = exp(A*X); 942 f3.print_symbolic(); 943 f3.print(); 944 } 945 946 TEST_CASE("testing function convexity"){ 947 var<> dp("dp",0.5,10.); 948 var<int> ip("ip",3,4); 949 auto expr = log(dp) + sqrt(ip); 950 expr.print_symbolic(); 951 CHECK(expr.is_concave()); 952 CHECK(expr._range->first==log(0.5)+(int)sqrt(3)); 953 CHECK(expr._range->second==log(10)+(int)sqrt(4)); 954 CHECK(expr.is_positive()); 955 var<> p("p",0,1); 956 var<> q("q",-0.5, 0.5); 957 auto cc = p*p + q*q; 958 cc.print_symbolic(); 959 CHECK(cc.is_convex()); 960 auto cc1 = cc * -1; 961 cc1.print_symbolic(); 962 CHECK(cc1.is_concave()); 963 cc1 += expr; 964 cc1.print_symbolic(); 965 CHECK(cc1.is_concave()); 966 param<int> aa("aa"); 967 aa = -1; 968 aa = -3; 969 auto ff = (aa)*p*p; 970 ff.print_symbolic(); 971 CHECK(ff.is_concave()); 972 ff *= aa; 973 ff.print_symbolic(); 974 CHECK(ff.is_convex()); 975 ff *= -1; 976 ff.print_symbolic(); 977 CHECK(ff.is_concave()); 978 ff *= aa; 979 ff.print_symbolic(); 980 CHECK(ff.is_convex()); 981 ff += aa*(ip + dp)*q*q; 982 CHECK(ff._range->first==-10.5); 983 CHECK(ff._range->second==37.5); 984 ff.print_symbolic(); 985 CHECK(!ff.is_convex()); 986 CHECK(!ff.is_concave()); 987 CHECK(ff.is_polynomial()); 988 CHECK(ff.get_nb_vars()==4); 989 param<> b("b"); 990 b = 1; 991 auto fn = p*p + q*q; 992 fn.print_symbolic(); 993 CHECK(fn.is_convex()); 994 fn += exp(p); 995 fn.print_symbolic(); 996 CHECK(fn.is_convex()); 997 auto f2 = pow(p,4); 998 CHECK(f2.is_convex()); 999 var<> x("x",1,2), y("y",2,4), z("z",2,4); 1000 x.in(R(1)); 1001 y.in(R(1)); 1002 z.in(R(1)); 1003 auto f3 = x*(0.01*(x-y)-1); 1004 f3.print(); 1005 auto f4 = x*(0.01*(x-y)-1) + z; 1006 f4.print(); 1007 CHECK(!f4.is_rotated_soc()); 1008 } 1009 1010 TEST_CASE("testing bounds copy"){ 1011 var<int> x("x",-3,3); 1012 var<> y("y", -1.2, 1.4); 1013 y.copy_bounds(x); 1014 y.print(); 1015 CHECK(y.get_lb(0)==-3); 1016 CHECK(y.get_ub(0)==3); 1017 } 1018 1019 TEST_CASE("testing bounds get"){ 1020 indices ids("index_set"); 1021 ids.add({"id1", "id2"}); 1022 param<> lb("lb"); 1023 lb.in(ids); 1024 lb.print(); 1025 lb("id1") = 1.5; 1026 lb("id2") = -231.5; 1027 param<> ub("ub"); 1028 ub.in(ids); 1029 ub.print(); 1030 ub("id1") = 2.5; 1031 ub("id2") = 231.5; 1032 var<> x("x",lb,ub); 1033 x.in(ids); 1034 CHECK(x.get_lb("id1")==1.5); 1035 CHECK(x.get_lb("id2")==-231.5); 1036 CHECK(x.get_ub(0)==2.5); 1037 CHECK(x.get_ub(1)==231.5); 1038 Constraint<> C("C"); 1039 C += x*x.get_lb(); 1040 C <= x.get_ub(); 1041 C.print(); 1042 x("id1").set_lb(2); 1043 x("id1").set_ub(2.1); 1044 x.set_lb("id2",-1); 1045 x.set_ub("id2",1); 1046 C.print();/* C has new params */ 1047 auto inst = 0, precision=10; 1048 auto str =C.to_str(inst,precision); 1049 CHECK(C.to_str(inst,precision)=="2x[id1] - 2.1"); 1050 } 1051 1052 TEST_CASE("testing quadratic function factorization"){ 1053 var<> x1("x1",0.5,10), x2("x2",-1,1); 1054 x1.in(R(4));x2.in(R(4)); 1055 auto f = 3*pow(x1,2) + 5*pow(x2,2) - 2*x1*x2; /* Convex since it can be factorized to (x1 - x2)^2 + 2x1^2 + 4x2^2 */ 1056 f.print_symbolic(); 1057 f.print(); 1058 CHECK(f.is_convex()); 1059 f = 4*x1*x2 - 2*pow(x1,2) - 3*pow(x2,2); /* Concave since it can be factorized to -(2-sqrt(2))x1^2 - (3-sqrt(2))x2^2 - (sqrt(2)x1 - sqrt(2)x2)^2*/ 1060 f.print_symbolic(); 1061 f.print(); 1062 CHECK(f.is_concave()); 1063 /* Checking convexity in model objective */ 1064 var<> x("x",-10,10); 1065 auto test = make_shared<Model<>>("test"); 1066 test->add(x.in(R(2))); 1067 x.initialize_uniform(); 1068 test->min(6*x(0)*x(0) + 4*x(1)*x(1) - 2.5*x(0)*x(1)); 1069 Constraint<> c1("c1"); 1070 c1 = x(0)*x(1); 1071 test->add(c1>=8); 1072 test->print(); 1073 solver<> s(test,ipopt); 1074 s.run(5, 1e-6); 1075 test->print_solution(); 1076 CHECK(std::abs(test->_obj->get_val()-58.3836718)<1e-6); 1077 } 1078 1079 TEST_CASE("testing constraints"){ 1080 var<> x1("x1", -1, 1), x2("x2", 0.1, 3), x3("x3",2,4); 1081 x1.in(R(2)); 1082 x2.in(R(2)); 1083 x3.in(R(2)); 1084 Constraint<> cstr("cstr"); 1085 cstr = x1 + exp(x2) - x3; 1086 CHECK(cstr.get_nb_vars()==3); 1087 CHECK(cstr.is_nonlinear()); 1088 CHECK(cstr.is_convex()); 1089 CHECK(cstr.get_dim()==2); 1090 CHECK(cstr._range->first==-1+exp(0.1)-4); 1091 CHECK(cstr._range->second==1+exp(3)-2); 1092 cstr.print_symbolic(); 1093 cstr.print(); 1094 auto dfdx2 = cstr.get_derivative(x2); 1095 dfdx2.print_symbolic(); 1096 param<int> a("a"); 1097 a = -1; 1098 a = -4; 1099 Constraint<> cstr1("cstr1"); 1100 cstr1 = x3 - sqrt(x2) + ReLU(x1); 1101 cstr1 >= a + 3; 1102 cstr1.print_symbolic(); 1103 cstr1.print(); 1104 CHECK(cstr1.is_concave()); 1105 CHECK(cstr1._range->first==2-sqrt(3)-3+1); 1106 CHECK(cstr1._range->second==4-sqrt(0.1)+2); 1107 } 1108 1109 TEST_CASE("testing soc/rotated soc constraints"){ 1110 var<> x1("x1", -1, 1), x2("x2", 0.1, 3), x3("x3",2,4); 1111 x1.in(R(2)); 1112 x2.in(R(2)); 1113 x3.in(R(2)); 1114 param<> a("a"); 1115 a = 1; 1116 a = 4; 1117 Constraint<> cstr("cstr"); 1118 cstr = a*x2*x3 - pow(x1,2); 1119 cstr >= 0; 1120 cstr.print_symbolic(); 1121 cstr.print(); 1122 CHECK(cstr.to_str()==" - x1² + (a)x2x3"); 1123 CHECK(cstr.get_nb_vars()==3); 1124 CHECK(cstr.is_quadratic()); 1125 CHECK(cstr.check_rotated_soc()); 1126 CHECK(cstr.get_dim()==2); 1127 auto dfdx2 = cstr.get_derivative(x2); 1128 dfdx2.print_symbolic(); 1129 CHECK(dfdx2.to_str()=="(a)x3"); 1130 CHECK(dfdx2.is_linear()); 1131 CHECK(dfdx2._range->first==2); 1132 CHECK(dfdx2._range->second==16); 1133 Constraint<> cstr2("cstr2"); 1134 cstr2 = pow(x2,2) + pow(a*x3,2) - pow(x1,2); 1135 cstr2 <= 0; 1136 cstr2.print_symbolic(); 1137 cstr2.print(); 1138 CHECK(cstr2.get_nb_vars()==3); 1139 CHECK(cstr2.is_quadratic()); 1140 CHECK(cstr2.check_soc()); 1141 CHECK(cstr2.get_dim()==2); 1142 CHECK(cstr2._range->first==0.01 + 2*2 - 1); 1143 CHECK(cstr2._range->second==3*3 + 16*4*4); 1144 } 1145 1146 TEST_CASE("testing nonlinear expressions"){ 1147 var<> x1("x1", -1, 1), x2("x2", 0.1, 3), x3("x3"); 1148 x1.in(R(1)); 1149 x2.in(R(1)); 1150 x3.in(R(1)); 1151 param<> a("a"), b("b"); 1152 a = 0.2;b = -0.6; 1153 param<> c("c"), d("d"); 1154 c = 0.1;d = -0.5; 1155 auto cstr = x1*exp(x2*x3); 1156 CHECK(cstr.get_nb_vars()==3); 1157 CHECK(cstr.is_nonlinear()); 1158 CHECK(cstr.get_dim()==1); 1159 cstr.print_symbolic(); 1160 cstr.print(); 1161 auto dfdx2 = cstr.get_derivative(x2); 1162 dfdx2.print_symbolic(); 1163 auto f = 2*cos(min(acos(a/b), asin(c/d)))*x1; 1164 f.print_symbolic(); 1165 CHECK(f.to_str()=="(2cos(min(acos(a/b), asin(c/d))))x1"); 1166 f.print(); 1167 } 1168 1169 TEST_CASE("testing monomials"){ 1170 var<> x1("x1"), x2("x2"), x3("x3"); 1171 auto cstr = x1*x2 + x2*x3 + x1*x3; 1172 auto monoms = cstr.get_monomials(5); /* Get all monomials up to degree 5 involving the variables in cstr */ 1173 CHECK(monoms.size()==21); 1174 } 1175 1176 TEST_CASE("testing compositions"){ 1177 auto v = build_compositions(4, 2); 1178 /* Should return 1179 3 1 1180 2 2 1181 1 3 1182 0 4 */ 1183 for (auto i = 0; i< v.size(); i++) { 1184 for (auto &row:v[i]) { 1185 cout << row << " "; 1186 } 1187 cout << endl; 1188 } 1189 } 1190 1191 TEST_CASE("testing simple model"){ 1192 Model<> M("MyModel"); 1193 var<> x("x", 1,2), y("y", 2, 4), z("z", -1, 2); 1194 M.add(x.in(R(2)),y.in(R(2)),z.in(R(2))); 1195 Constraint<> cstr1("cstr1"); 1196 cstr1 = pow(x,2) + pow(y,2) - pow(z,2); 1197 M.add(cstr1 <= 0); 1198 Constraint<> cstr2("cstr2"); 1199 cstr2 = 1; 1200 M.add(cstr2 >= 2); 1201 M.max(sum(x)); 1202 M.print(); 1203 M.round_solution(); 1204 CHECK(M.is_convex()); 1205 CHECK(M.get_nb_vars()==6); 1206 CHECK(M.get_nb_cons()==2); 1207 } 1208 1209 TEST_CASE("testing nonlinear Model"){ 1210 Model<> M("MyModel2"); 1211 param<> x_lb("x_lb"), x_ub("x_ub"), y_lb("y_lb"), y_ub("y_ub"); 1212 x_lb.in(R(4));x_ub.in(R(4));y_lb.in(R(4));y_ub.in(R(4)); 1213 x_lb = -1; x_ub = 4; y_lb = -5; y_ub = 3; 1214 x_lb(2) = -12; 1215 var<> x("x", x_lb,x_ub), y("y", y_lb, y_ub), z("z", -1, 2); 1216 M.add(x.in(R(4)),y.in(R(4)),z.in(R(4))); 1217 CHECK(x.get_lb(2)==-12); 1218 Constraint<> cstr1("cstr1"); 1219 cstr1 = pow(x,2) + pow(y,2) - pow(z,2); 1220 M.add(cstr1 <= 0); 1221 param<> a("a"); 1222 a = 2;a = 3;a = 4;a = 5; 1223 Constraint<> cstr2("cstr2"); 1224 cstr2 = a*x*y*cos(x-z); 1225 M.add(cstr2 == 2); 1226 CHECK(cstr2.to_str()==" - 2 + (a)xy * cos(x - z)"); 1227 M.max(sum(x)); 1228 M.print_symbolic(); 1229 M.print(); 1230 CHECK(!M.is_convex()); 1231 CHECK(M.get_nb_vars()==12); 1232 CHECK(M.get_nb_cons()==8); 1233 } 1234 1235 TEST_CASE("testing complex constraint expansion"){ 1236 Model<> M("test"); 1237 auto ids = indices("ids"); 1238 ids.add({"id1", "id2"}); 1239 var<> x("x",-1,1), y("y",-2,2), u1("u1",0,1), v1("v1",-3,1), u2("u2",2,3), v2("v2",0.5,1.4); 1240 M.add(x.in(ids),y.in(ids),u1.in(ids),v1.in(ids),u2.in(ids),v2.in(ids)); 1241 var<Cpx> w1("w1"), w2("w2"), z("z"); 1242 z.real_imag(x, y); 1243 w1.real_imag(u1, v1); 1244 w2.real_imag(u2, v2); 1245 1246 param<> pr1("pr1"), pi1("pi1"),pr2("pr2"); 1247 pr1 = {1,2};pi1 = {0,-1};pr2 = {-2,2}; 1248 param<Cpx> p1("p1"); 1249 p1.real_imag(pr1, pi1); 1250 1251 Constraint<Cpx> C_lin1("C_lin1"); 1252 C_lin1 = p1*z; 1253 M.add(C_lin1.in(ids)==0); 1254 1255 param<Cpx> p2("p2"); 1256 p2.set_real(pr2);/* zero imaginary */ 1257 Constraint<Cpx> C_lin2("C_lin2"); 1258 C_lin2 = p2*z; 1259 M.add(C_lin2.in(ids)==0); 1260 1261 Constraint<Cpx> C_quad("C_quad"); 1262 C_quad = z*w1; 1263 M.add(C_quad.in(ids)==0); 1264 1265 Constraint<Cpx> C_norm("C_norm"); 1266 C_norm = z*conj(z); 1267 M.add(C_norm.in(ids)==0); 1268 1269 Constraint<Cpx> C_pol("C_pol"); 1270 C_pol = z*w1*w2; 1271 M.add(C_pol.in(ids)==0); 1272 M.print(); 1273 } 1274 1275 TEST_CASE("testing multithread solve"){ 1276 string fname = string(prj_dir)+"/data_sets/Power/nesta_case5_pjm.m"; 1277 unsigned nb_threads = 2; 1278 double tol = 1e-6; 1279 string mehrotra = "no", log_level="0"; 1280 PowerNet grid1,grid2; 1281 grid1.readgrid(fname); 1282 auto ACOPF1 = build_ACOPF(grid1,ACRECT); 1283 fname = string(prj_dir)+"/data_sets/Power/nesta_case9_bgm__nco.m"; 1284 grid2.readgrid(fname); 1285 auto ACOPF2 = build_ACOPF(grid2,ACPOL); 1286 auto models = {ACOPF1, ACOPF2}; 1287 /* run in parallel */ 1288 // solver<> OPF1(ACOPF1,ipopt); 1289 // solver<> OPF2(ACOPF2,ipopt); 1290 // OPF1.run(0,tol); 1291 // OPF2.run(0,tol); 1292 run_parallel(models, ipopt, tol = 1e-6, nb_threads=2); 1293 CHECK(std::abs(ACOPF1->get_obj_val()-17551.89092)<1e-2); 1294 CHECK(std::abs(ACOPF2->get_obj_val()-3087.83977)<1e-2); 1295 CHECK(ACOPF1->is_feasible(tol)); 1296 ACOPF1->print_solution(); 1297 auto Mc = ACOPF1->build_McCormick(); 1298 auto clone = grid1.clone(); 1299 CHECK(clone->arcs.size()==grid1.arcs.size()); 1300 clone->remove_arc(clone->arcs.at(0)); 1301 CHECK(clone->arcs.size()==grid1.arcs.size()-1); 1302 } 1303 1304 TEST_CASE("testing SOCOPF anc Cycle basis computation"){ 1305 string fname = string(prj_dir)+"/data_sets/Power/nesta_case5_pjm.m"; 1306 // string fname = "/Users/hlh/Dropbox/Work/Dev/pglib-opf-18.08/pglib_opf_case2383wp_k.m"; 1307 int output = 0; 1308 double tol = 1e-6; 1309 string mehrotra = "no", log_level="0"; 1310 PowerNet grid; 1311 grid.readgrid(fname); 1312 auto SOCOPF = grid.build_SCOPF(); 1313 // SOCOPF->print(); 1314 CHECK(SOCOPF->_type==quad_m); 1315 solver<> OPF(SOCOPF,ipopt); 1316 auto time_start = get_wall_time(); 1317 OPF.run(output=5, tol=1e-6); 1318 auto time_end = get_wall_time(); 1319 DebugOn("Total wall time = " << time_end - time_start << " secs" << endl); 1320 CHECK(std::abs(SOCOPF->_obj->get_val()-14999.715)<1e-2); 1321 auto cycle_basis = grid.get_cycle_basis();/** Compute cycle basis */ 1322 for (auto cycle: cycle_basis) { 1323 cycle->print(); 1324 } 1325 } 1326 1327 TEST_CASE("testing SDPOPF"){ 1328 string fname = string(prj_dir)+"/data_sets/Power/pglib_opf_case24_ieee_rts__api.m"; 1329 int output = 0; 1330 double tol = 1e-6; 1331 string mehrotra = "no", log_level="0"; 1332 PowerNet grid; 1333 grid.readgrid(fname); 1334 bool add_current_RLT = true; 1335 auto OPF = build_ACOPF(grid, ACRECT); 1336 auto SDP = build_SDPOPF(grid, add_current_RLT); 1337 solver<> UB(OPF,ipopt); 1338 solver<> LB(SDP,ipopt); 1339 auto time_start = get_wall_time(); 1340 UB.run(output=5, tol=1e-6); 1341 LB.run(output=5, tol=1e-6); 1342 auto final_gap = 100*(OPF->get_obj_val() - SDP->get_obj_val())/std::abs(OPF->get_obj_val()); 1343 auto time_end = get_wall_time(); 1344 DebugOn("Total wall time = " << time_end - time_start << " secs" << endl); 1345 DebugOn("SDP Gap = " << final_gap << "%" << endl); 1346 CHECK(final_gap<1.2); 1347 } 1348 1349 TEST_CASE("Bug in Cplex MIQCP presolve"){ 1350 /* This test identifies a bug in Cplex's 12.8 MIQCP presolve. Turn Cplex = true and relax = false to get an infeasible status, with relax=true, Cplex returns feasible. Also turning off presolve returns feasible. */ 1351 1352 /* Start indexing from 1 */ 1353 indices R1 = range(1,1); 1354 indices R2 = range(1,2); 1355 indices R4 = range(1,4); 1356 indices R10 = range(1,10); 1357 1358 Model<> m; 1359 1360 /* Bounds */ 1361 param<> x_lb("x_lb"); 1362 param<> x_ub("x_ub"); 1363 x_lb = {-2,-2,-1,35,-4,-3,0,0,-4,0}; 1364 x_ub = {2,2,1,37,-2,5,25,4,4,4}; 1365 1366 var<> x("x",x_lb,x_ub); 1367 var<int> A1("A1",0,1), A2("A2",0,1), A6("A6",0,1); 1368 m.add(A1.in(R1), A2.in(R1), A6.in(R1)); 1369 var<> L7("L7",0,1), L8("L8",0,1), L9("L9",0,1), L10("L10",0,1); 1370 m.add(x.in(R10), L7.in(R2), L8.in(R2), L9.in(R4), L10.in(R2)); 1371 1372 1373 m.min(x[3]+x[4]+x[5]); 1374 1375 Constraint<> c1("c1"); 1376 c1 = x[3] - x[7]; 1377 m.add(c1==0); 1378 Constraint<> c2("c2"); 1379 c2 = x[4] + 14* x[1] - 3*x[8] + 14*x[2] - 6*x[9] - 3*x[10] - 19; 1380 m.add(c2==0); 1381 Constraint<> c3("c3"); 1382 c3 = x[5] + 32* x[1] - 12*x[8] - 48*x[2] + 36*x[9] - 27*x[10] - 18; 1383 m.add(c3==0); 1384 Constraint<> c4("c4"); 1385 c4 = x[6] - x[2] - x[1] - 1; 1386 m.add(c4==0); 1387 Constraint<> c5("c5"); 1388 c5 = A1[1] - 1; 1389 m.add(c5==0); 1390 Constraint<> c6("c6"); 1391 c6 = x[1] + 2*A1[1]; 1392 m.add(c6>=0); 1393 Constraint<> c7("c7"); 1394 c7 = x[1] - 2*A1[1]; 1395 m.add(c7<=0); 1396 Constraint<> c8("c8"); 1397 c8 = L8[1] + L8[2]; 1398 m.add(c8==1); 1399 Constraint<> c9("c9"); 1400 c9 = x[8] - 4*L8[1] - 4* L8[2]; 1401 m.add(c9<=0); 1402 Constraint<> c10("c10"); 1403 c10 = L8[1] - A1[1]; 1404 m.add(c10<=0); 1405 Constraint<> c11("c11"); 1406 c11 = L8[2] - A1[1]; 1407 m.add(c11<=0); 1408 Constraint<> c12("c12"); 1409 c12 = x[1] + 2* A1[1]; 1410 m.add(c12>=0); 1411 Constraint<> c13("c13"); 1412 c13 = x[1] - 2* A1[1]; 1413 m.add(c13<=0); 1414 Constraint<> c14("c14"); 1415 c14 = x[1] + 2* L8[1] - 2* L8[2]; 1416 m.add(c14==0); 1417 Constraint<> c15("c15"); 1418 c15 = A2[1] - 1; 1419 m.add(c15==0); 1420 Constraint<> c16("c16"); 1421 c16 = x[2] + 2* A2[1]; 1422 m.add(c16>=0); 1423 Constraint<> c17("c17"); 1424 c17 = x[2] - 2* A2[1]; 1425 m.add(c17<=0); 1426 Constraint<> c18("c18"); 1427 c18 = L9[1] + L9[2] + L9[3] + L9[4] - 1; 1428 m.add(c18==0); 1429 Constraint<> c19("c19"); 1430 c19 = x[9] - 4* L9[1] + 4* L9[2] + 4* L9[3] - 4* L9[4]; 1431 m.add(c19==0); 1432 Constraint<> c20("c20"); 1433 c20 = L9[1] + L9[3] - A2[1]; 1434 m.add(c20<=0); 1435 Constraint<> c21("c21"); 1436 c21 = L9[2] + L9[4] - A2[1]; 1437 m.add(c21<=0); 1438 Constraint<> c22("c22"); 1439 c22 = x[2] + 2* L9[1] + 2* L9[3] - 2* L9[2] - 2* L9[4]; 1440 m.add(c22==0); 1441 Constraint<> c23("c23"); 1442 c23 = L9[1] + L9[2] - A1[1]; 1443 m.add(c23<=0); 1444 Constraint<> c24("c24"); 1445 c24 = L9[3] + L9[4] - A1[1]; 1446 m.add(c24<=0); 1447 Constraint<> c25("c25"); 1448 c25 = x[1] + 2* L9[1] + 2* L9[2] - 2* L9[3] - 2* L9[4]; 1449 m.add(c25==0); 1450 Constraint<> c26("c26"); 1451 c26 = L10[1] + L10[2]; 1452 m.add(c26==1); 1453 Constraint<> c27("c27"); 1454 c27 = x[10] - 4* L10[1] - 4* L10[2]; 1455 m.add(c27<=0); 1456 Constraint<> c28("c28"); 1457 c28 = x[2] + 2* L10[1] - 2* L10[2]; 1458 m.add(c28==0); 1459 Constraint<> c29("c29"); 1460 c29 = A6[1]; 1461 m.add(c29==1); 1462 Constraint<> c30("c30"); 1463 c30 = L7[1] + L7[2]; 1464 m.add(c30==1); 1465 Constraint<> c31("c31"); 1466 c31 = x[7] - 9* L7[1] - 25* L7[2]; 1467 m.add(c31<=0); 1468 Constraint<> c32("c32"); 1469 c32 = x[6] + 3* L7[1] - 5* L7[2]; 1470 m.add(c32==0); 1471 Constraint<> c33("c33"); 1472 c33 = -1*pow(x[1],2) + x[8]; 1473 m.add(c33>=0); 1474 Constraint<> c34("c34"); 1475 c34 = -1*pow(x[2],2) + x[10]; 1476 m.add(c34>=0); 1477 Constraint<> c35("c35"); 1478 c35 = -1*pow(x[6],2) + x[7]; 1479 m.add(c35>=0); 1480 m.print(); 1481 bool use_cplex = false, relax = true; 1482 if(use_cplex){ 1483 solver<> s(m,cplex); 1484 s.run(relax); 1485 } 1486 else { 1487 solver<> s(m,ipopt); 1488 s.run(5,1e-6); 1489 } 1490 CHECK(std::abs(m.get_obj_val()-31.00313)<1e-3); 1491 m.print_solution(); 1492 } 1493 1494 TEST_CASE("Alpine issue") { 1495 Model<> m; 1496 var<> x1("x1", -2, 2), x2("x2", -2, 2); 1497 var<> y1("y1"), y2("y2"), y3("y3"), y4("y4"); 1498 m.add(x1.in(R(1)), x2.in(R(1)),y1.in(R(1)), y2.in(R(1)),y3.in(R(1)),y4.in(R(1))); 1499 Constraint<> c1("c1"); 1500 c1 = pow((x1 +x2 +1),2); 1501 CHECK(c1.is_convex()); 1502 y1.set_lb(c1._range->first); 1503 y1.set_ub(c1._range->second); 1504 c1 -= y1; 1505 m.add(c1==0); 1506 cout << "y1 lower bound = " << y1._range->first << endl; 1507 cout << "y1 upper bound = " << y1._range->second << endl; 1508 Constraint<> c2("c2"); 1509 // c2 = 19 - 14*x1 + 3*pow(x1,2) - 14*x2 + 6*x1*x2 + 3*pow(x2,2); 1510 c2 = 19 - 14*x1 - 14*x2 + pow(sqrt(3)*x1+sqrt(3)*x2,2); 1511 y2.set_lb(c2._range->first); 1512 y2.set_ub(c2._range->second); 1513 cout << "y2 lower bound = " << y2._range->first << endl; 1514 cout << "y2 upper bound = " << y2._range->second << endl; 1515 c2 -= y2; 1516 CHECK(c2.is_convex()); 1517 m.add(c2==0); 1518 Constraint<> c3("c3"); 1519 c3 = pow(2*x1 - 3*x2,2); 1520 y3.set_lb(c3._range->first); 1521 y3.set_ub(c3._range->second); 1522 c3 -= y3; 1523 CHECK(c3.is_convex()); 1524 m.add(c3==0); 1525 cout << "y3 lower bound = " << y3._range->first << endl; 1526 cout << "y3 upper bound = " << y3._range->second << endl; 1527 Constraint<> c4("c4"); 1528 // c4 = 18 - 32*x1 + 12*pow(x1,2) + 48*x2 - 36*x1*x2 + 27*pow(x2,2); 1529 c4 = 18 - 32*x1 - 6*pow(x1,2) + 48*x2 + pow(sqrt(18)*x1 - sqrt(18)*x2,2) + 9*pow(x2,2); 1530 y4.set_lb(c4._range->first); 1531 y4.set_ub(c4._range->second); 1532 c4 -= y4; 1533 CHECK(!c4.is_convex()); 1534 m.add(c4==0); 1535 cout << "y4 lower bound = " << y4._range->first << endl; 1536 cout << "y4 upper bound = " << y4._range->second << endl; 1537 1538 auto obj = 30 + y3*y4 +30*y1*y2 + y1*y2*y3*y4; 1539 cout << "Objective lower bound = " << obj._range->first << endl; 1540 cout << "Objective upper bound = " << obj._range->second << endl; 1541 m.min(obj); 1542 m.print_symbolic(); 1543 // m.initialize_uniform(); 1544 solver<> s(m,ipopt); 1545 // s.run(5,1e-6,1000000000); 1546 // m.print_solution(); 1547 Model<> relax("Lifted Relaxation"); 1548 relax.add(x1.in(R(1)), x2.in(R(1)),y1.in(R(1)), y2.in(R(1)),y3.in(R(1)),y4.in(R(1))); 1549 var<> y1234("y1234"), y12("y12"), y34("y34"), x12("(x12"), x11("x11"), x22("x22"); 1550 relax.add(x11.in(R(1)), x22.in(R(1)),x12.in(R(1)), y12.in(R(1)),y34.in(R(1)),y1234.in(R(1))); 1551 y1234.set_lb((y1*y2*y3*y4)._range->first); 1552 y1234.set_ub((y1*y2*y3*y4)._range->second); 1553 y12.set_lb((y1*y2)._range->first); 1554 y12.set_ub((y1*y2)._range->second); 1555 y34.set_lb((y3*y4)._range->first); 1556 y34.set_ub((y3*y4)._range->second); 1557 x12.set_lb((x1*x2)._range->first); 1558 x12.set_ub((x1*x2)._range->second); 1559 x11.set_lb((x1*x1)._range->first); 1560 x11.set_ub((x1*x1)._range->second); 1561 x22.set_lb((x1*x1)._range->first); 1562 x22.set_ub((x2*x2)._range->second); 1563 // func<> relax_obj = 30 + y34 +30*y12 + y1234; 1564 func<> relax_obj = y1234; 1565 relax.max(relax_obj); 1566 relax.add(c1<=0); 1567 relax.add(c2<=0); 1568 relax.add(c3<=0); 1569 1570 Constraint<> obj_ub("obj_ub"); 1571 1572 obj_ub = 30 + y34 +30*y12 + y1234; 1573 relax.add(obj_ub<=3); 1574 1575 Constraint<> c1_relax("c1_relax"); 1576 c1_relax = x11 + 2*x12 + x22 + 2*x1 + 2*x2 - y1 + 1; 1577 relax.add(c1_relax==0); 1578 Constraint<> c2_relax("c2_relax"); 1579 c2_relax = 3*x11 + 6*x12 + 3*x22 - 14*x1 - 14*x2 - y2 + 19; 1580 relax.add(c2_relax==0); 1581 Constraint<> c3_relax("c3_relax"); 1582 c3_relax = 4*x11 - 12*x12 + 9*x22 - y3; 1583 relax.add(c3_relax==0); 1584 Constraint<> c4_relax("c4_relax"); 1585 c4_relax = 12*x11 -36*x12 + 27*x22 - 32*x1 + 48*x2 - y4 + 18; 1586 relax.add(c4_relax==0); 1587 1588 Constraint<> soc1("soc1"); 1589 soc1 = x11 - pow(x1,2); 1590 relax.add(soc1 >= 0); 1591 Constraint<> soc2("soc2"); 1592 soc2 = x22 - pow(x2,2); 1593 relax.add(soc2 >= 0); 1594 1595 Constraint<> rot_soc("rot_soc"); 1596 rot_soc = x11*x22 - pow(x12,2); 1597 relax.add(rot_soc >= 0); 1598 CHECK(rot_soc.is_rotated_soc()); 1599 1600 relax.add_McCormick("x12", x12, x1, x2); 1601 relax.add_McCormick("y12", y12, y1, y2); 1602 relax.add_McCormick("y34", y34, y3, y4); 1603 relax.add_McCormick("y1234", y1234, y12, y34); 1604 CHECK(relax.is_convex()); 1605 relax.print_symbolic(); 1606 solver<> s2(relax,ipopt); 1607 s2.run(); 1608 relax.print_solution(); 1609 } 1610 1611 TEST_CASE("testing absolute value function") { 1612 1613 Model<> M("Test"); 1614 var<> x("x", -6, 2); 1615 var<> y("y", -3, 2); 1616 var<> obj("obj", pos_); 1617 M.add(x.in(R(1)),y.in(R(1)),obj.in(R(1))); 1618 1619 M.min(abs(2*x)*y + 36); 1620 1621 M.print(); 1622 solver<> NLP(M,ipopt); 1623 int output; 1624 double tol; 1625 string lin_solver; 1626 NLP.run(output=5,tol=1e-6); 1627 M.print_solution(); 1628 M.print_symbolic(); 1629 CHECK(M.get_obj_val()<=tol); 1630 } 1631 1632 TEST_CASE("testing min/max functions") { 1633 Model<> M("Test"); 1634 var<> x("x", -6, 3); 1635 var<> y("y", -3, 2); 1636 var<> obj("obj", -10,10); 1637 M.add(x.in(R(1)),y.in(R(1)),obj.in(R(1))); 1638 1639 M.min(min(x,3*y) + 6 - max(2*x,y)); 1640 1641 M.print(); 1642 solver<> NLP(M,ipopt); 1643 int output; 1644 double tol; 1645 string lin_solver; 1646 NLP.set_option("bound_relax_factor", 1e-10); 1647 NLP.run(output=5,tol=1e-6); 1648 M.print_solution(); 1649 M.print_symbolic(); 1650 CHECK(std::abs(M.get_obj_val()-(-9))<tol); 1651 } 1652 1653 1654 TEST_CASE("testing normal distributions") { 1655 param<> A("A"); 1656 int n = 200, p = 50; 1657 A.in(R(n,p)); 1658 double mean = 0, dev = 1; 1659 A.initialize_normal(mean, dev); 1660 1661 Model<> M("Normal"); 1662 var<> X("X", -1, 1), Xp("Xp", pos_); 1663 M.add(X.in(R(p)), Xp.in(R(p))); 1664 1665 var<> s("slack"); 1666 M.add(s.in(R(n))); 1667 1668 // M.min(sum(Xp) + 1e3*product(1,pow(s,2))); 1669 M.min(sum(Xp) + 1e3*norm2(s)); 1670 1671 Constraint<> Abs_p("x_abs_p"); 1672 Abs_p = X - Xp; 1673 M.add(Abs_p <= 0); 1674 1675 Constraint<> Abs_n("x_abs_n"); 1676 Abs_n = Xp + X; 1677 M.add(Abs_n >= 0); 1678 1679 Constraint<> Norm2("norm2"); 1680 // Norm2 = product(1,pow(X,2)); 1681 Norm2 = norm2(X); 1682 M.add(Norm2==1); 1683 1684 Constraint<> Lin("lin"); 1685 Lin = product(A,X); 1686 M.add(Lin==s); 1687 1688 M.print_symbolic(); 1689 } 1690 1691 TEST_CASE("testing set union unindexed") { 1692 indices ids1("index_set1"); 1693 ids1 = indices(range(1,2), range(2,4)); 1694 indices ids2("index_set2"); 1695 ids2 = indices(range(1,5), range(2,4)); 1696 auto union_set = union_ids(ids1, ids2); 1697 union_set.print(); 1698 CHECK(union_set.size()==15); 1699 indices ids3("index_set3"); 1700 ids3 = indices(range(1,5)); 1701 REQUIRE_THROWS_AS(union_ids(ids1,ids3), invalid_argument); 1702 } 1703 1704 TEST_CASE("testing set union indexed") { 1705 DebugOn("testing set union indexed" << endl); 1706 indices ids1("index_set1"); 1707 ids1.add("1","2","3","4","5","6","7","8","9"); 1708 var<> v1("v1"); 1709 v1.in(ids1); 1710 1711 indices ids2("index_set2"); 1712 ids2.add("5,4", "3,5", "6,7", "7,8", "2,6", "8,9"); 1713 var<> v2("v2"); 1714 v2 = v1.from(ids2); 1715 1716 var<> v3("v3"); 1717 v3 = v1.to(ids2); 1718 1719 auto union_set = union_ids(*v2._indices, *v3._indices); //in this case, the keys are same, so the union should check not only keys but also ._ids in the individual index sets and add accordingly 1720 CHECK(union_set.size()==8); //in the _ids, the function should work in a way that union_set._ids = [4,2,5,6,1,7,3,8] 1721 REQUIRE_THROWS_AS(union_ids(ids1,ids2), invalid_argument); 1722 } 1723 1724 TEST_CASE("testing from_ith() function") { 1725 indices ids("index_set"); 1726 ids = indices(range(1,3),range(9,10), range(2,4)); 1727 param<> dp("dp"); 1728 dp.in(range(2,4)); 1729 dp.print(); 1730 dp("2") = 1.5; 1731 dp("4") = -231.5; 1732 dp.print(); 1733 CHECK(dp.eval("2")==1.5); 1734 CHECK(dp.eval("4")==-231.5); 1735 CHECK(dp._range->first==-231.5); 1736 CHECK(dp._range->second==1.5); 1737 REQUIRE_THROWS_AS(dp("unexisting_key").eval(), invalid_argument); 1738 ids.print(); 1739 auto ndp = dp.from_ith(2,ids); 1740 ndp.print(); 1741 CHECK(ndp.get_dim()==ids.size()); 1742 indices ids2("index_set2"); 1743 ids2 = indices(range(1,3),range(9,10), range(2,4)); 1744 var<> dv("dv"); 1745 dv.in(range(9,10)); 1746 int precision = 5; 1747 dv.print_vals(precision=5); 1748 auto ndv = dv.from_ith(1,ids2); 1749 ndv.print_vals(precision=5); 1750 var<> dv2("dv2"); 1751 dv2.in(range(1,3)); 1752 auto ndv2 = dv2.from_ith(0,ids2); 1753 ndv2.print_vals(precision=5); 1754 } 1755 1756 1757 TEST_CASE("testing in_ignore_ith() function") { 1758 indices ids("index_set"); 1759 ids = indices(range(1,3),range(9,10), range(2,4)); 1760 param<> dp("dp"); 1761 dp.in(range(1,3),range(2,4)); 1762 dp("1,2") = 1.5; 1763 dp("3,4") = -231.5; 1764 dp.print(); 1765 auto ndp = dp.in_ignore_ith(1, 1, ids); 1766 ndp.print(); 1767 CHECK(ndp.get_dim()==ids.size()); 1768 indices idsv1("index_setv1"); 1769 idsv1.add("id1", "id11", "id111"); 1770 indices idsv2("index_setv2"); 1771 idsv2.add("id2", "id22", "id222"); 1772 var<> dv("dv"); 1773 dv.in(idsv1,idsv2); 1774 int precision = 5; 1775 dv.print_vals(precision=5); 1776 indices ids_all("index_set_all"); 1777 ids_all.add("id1,id3,id2", "id11,id33,id22", "id111,id333,id222"); 1778 auto ndv = dv.in_ignore_ith(1, 1, ids_all); 1779 ndv.print_vals(precision=5); 1780 } 1781 1782 TEST_CASE("testing get_matrix()") { 1783 auto ids = indices(range(1,3),range(8,12)); 1784 var<> dv("dv"); 1785 dv = dv.in(ids); 1786 dv.print_vals(4); 1787 1788 Constraint<> Sum1("Sum1"); 1789 Sum1 = sum(dv.in_matrix(0,1)); // when we say this it should start from entry 0 and get the first indices range(1,3) in the second dimension 1790 Sum1.print(); 1791 CHECK(Sum1.get_nb_instances() == 5); // so when we actually sum, the number of instances should be equal to range(8,12) = 5 1792 1793 Constraint<> Sum2("Sum2"); 1794 Sum2 = sum(dv.in_matrix(1,1)); 1795 Sum2.print(); 1796 CHECK(Sum2.get_nb_instances() == 3); 1797 1798 auto ids1 = indices(range(1,3),range(4,7),range(8,12)); 1799 var<> dv1("dv1"); 1800 dv1 = dv1.in(ids1); 1801 dv1.print_vals(4); 1802 1803 auto dv4 = dv1.in_matrix(1,1); //this case is wrong as well, it was working before because all the entries had only 3 elements. 1804 Constraint<> Sum3("Sum3"); 1805 Sum3 = sum(dv4); 1806 Sum3.print(); 1807 CHECK(Sum3.get_nb_instances() == 15); 1808 } 1809 1810 TEST_CASE("testing sum_ith()") { 1811 indices ids1("index_set1"); 1812 ids1.add("5,4,1", "5,2,1", "7,8,4", "5,5,1", "7,6,4", "7,9,4"); 1813 var<> v1("v1"); 1814 v1.in(ids1); 1815 Constraint<> Sum1("Sum1"); 1816 Sum1 = sum_ith(v1, 1, 1); 1817 stringstream buffer; 1818 auto console = cout.rdbuf(buffer.rdbuf()); 1819 Sum1.print(); 1820 cout.rdbuf(console); 1821 CHECK(buffer.str()==" Sum1 (Linear) : \nSum1[0]: v1[5,4,1] + v1[5,2,1] + v1[5,5,1] <= 0;\nSum1[1]: v1[7,8,4] + v1[7,6,4] + v1[7,9,4] <= 0;\n"); 1822 CHECK(Sum1.get_nb_instances() == 2); 1823 indices ids2("index_set2"); 1824 ids2.add("4,1", "2,1", "8,4"); 1825 indices ids3("index_set3"); 1826 ids3.add("4,1", "2,1"); 1827 param<> p2("p2"); 1828 p2.in(ids2); 1829 var<> v2("p2"); 1830 v2.in(ids2); 1831 p2("4,1") = -3.4; 1832 p2("8,4") = 1.5; 1833 p2("2,1") = -6; 1834 Constraint<> Sum2("Sum2"); 1835 Sum2 = sum(v2,ids3) + sum(p2,ids3); 1836 Sum2.print(); 1837 } 1838 1839 TEST_CASE("sum over outgoing") { 1840 /* Define indices */ 1841 indices arcs("arcs"); 1842 arcs.add("a1,1,2", "a2,1,3", "a3,1,4", "a4,3,4", "a5,2,4", "a6,1,5", "a7,1,5"); 1843 indices nodes("nodes"); 1844 nodes.add("1", "2", "3", "4", "5"); 1845 /** Declare model,vars and constraints */ 1846 Model<> M("Test"); 1847 var<> v1("flux", 0, 1); 1848 M.add(v1.in(arcs)); 1849 Constraint<> Sum0("Sum0"); 1850 Sum0 = v1.sum_out(nodes) + v1.sum_in(nodes); 1851 M.add(Sum0.in(nodes) == 0); 1852 M.print(); 1853 CHECK(Sum0.get_nb_instances() == nodes.size()); 1854 } 1855 1856 TEST_CASE("testing sum_ith() func<> version"){ 1857 indices ids1("index set1"""); 1858 ids1 = indices(range(1,3),range(1,4),range(1,6)); 1859 1860 indices ids2("index set2"); 1861 ids2 = indices(range(1,3),range(1,4),range(1,5),range(1,6)); 1862 1863 param<> p1("p1"); 1864 p1.in(ids1); 1865 size_t pos = 0; 1866 for(auto i = 1; i<= 3; i++){ 1867 for(auto j = 1; j<= 4; j++){ 1868 for(auto k = 1; k<= 6; k++){ 1869 p1.set_val(pos, 100*i+10*j+k); 1870 pos++; 1871 } 1872 } 1873 } 1874 1875 var<> v2("v2"); 1876 v2.in(ids2); 1877 auto pp1 = p1.in_ignore_ith(2,1,ids2); 1878 pp1.print_vals(4); 1879 Constraint<> Sum1("Sum1"); 1880 Sum1 = sum_ith(v2,1,2); 1881 CHECK(Sum1.get_nb_instances() == 3*6); 1882 1883 } 1884 1885 TEST_CASE("testing get_interaction_grap()") { 1886 indices buses("buses"); 1887 buses.insert("1", "2", "3", "4"); 1888 indices node_pairs("bpairs"); 1889 node_pairs.insert("1,2", "1,3", "3,4", "4,1"); 1890 1891 Model<> Mtest("Mtest"); 1892 var<> R_Wij("R_Wij", -1, 1); 1893 var<> Im_Wij("Im_Wij", -1, 1); 1894 var<> Wii("Wii", 0.8, 1.21); 1895 Mtest.add(R_Wij.in(node_pairs), Im_Wij.in(node_pairs), Wii.in(buses)); 1896 Constraint<> SOC("SOC"); 1897 SOC = pow(R_Wij, 2) + pow(Im_Wij, 2) - Wii.from(node_pairs)*Wii.to(node_pairs); 1898 Mtest.add(SOC.in(node_pairs) <= 0); 1899 Mtest.print(); 1900 auto g = Mtest.get_interaction_graph(); 1901 g.print(); 1902 CHECK(g.nodes.size()==12); 1903 CHECK(g.arcs.size()==4); 1904 } 1905 1906 1907 TEST_CASE("testing constraint delete") { 1908 indices buses("buses"); 1909 buses.insert("1", "2", "3", "4"); 1910 indices node_pairs("bpairs"); 1911 node_pairs.insert("1,2", "1,3", "3,4", "4,1"); 1912 1913 Model<> Mtest("Mtest"); 1914 var<> R_Wij("R_Wij", -1, 1); 1915 /* Imaginary part of Wij = ViVj */ 1916 var<> Im_Wij("Im_Wij", -1, 1); 1917 var<> Wii("Wii", 0.8, 1.21); 1918 Mtest.add(R_Wij.in(node_pairs), Im_Wij.in(node_pairs), Wii.in(buses)); 1919 Constraint<> SOC("SOC"); 1920 SOC = pow(R_Wij, 2) + pow(Im_Wij, 2) - Wii.from(node_pairs); 1921 Mtest.add(SOC.in(node_pairs)); 1922 1923 Constraint<> PAD("PAD"); 1924 PAD = 2*R_Wij - Im_Wij; 1925 Mtest.add(PAD.in(node_pairs)<=2); 1926 1927 Mtest.print(); 1928 CHECK(Mtest.get_nb_cons() == 8); 1929 Mtest.remove("SOC"); 1930 Mtest.print(); 1931 CHECK(Mtest.is_linear()); 1932 CHECK(Mtest.get_nb_cons() == 4); 1933 } 1934 1935 1936 1937 1938 1939 1940 TEST_CASE("Few Degrees Of Freedom") { 1941 auto m = Model<>("test"); 1942 /* ----- Indices ----- */ 1943 auto pair_ids = indices("pair_ids"); 1944 pair_ids.add("40,41", "41,42", "40,43", "41,43"); 1945 1946 1947 /* ----- Variables ----- */ 1948 auto x = var<>("x", -1, 1); 1949 m.add(x.in(range(40,43))); 1950 1951 /* ----- Objective ----- */ 1952 m.min(x[40] + x[41] - x[42] - x[43]); 1953 1954 /* ----- Constraints ----- */ 1955 auto Linear = Constraint<>("Linear"); 1956 Linear = x[40] - x[41]; 1957 m.add(Linear==0.5); 1958 1959 auto Bilinear = Constraint<>("Bilinear"); 1960 Bilinear = x.from(pair_ids)*x.to(pair_ids); 1961 m.add(Bilinear.in(pair_ids)==0); 1962 1963 m.print(); 1964 /* ----- Solver ----- */ 1965 auto s = solver<>(m,ipopt); 1966 s.run(5,1e-6); 1967 m.print_solution(); 1968 CHECK(abs(m.get_obj_val()-(-0.5))<1e-4); 1969 } 1970 1971 1972 TEST_CASE("Paths") { 1973 indices buses("buses"); 1974 buses.insert("1", "2", "3", "4"); 1975 indices node_pairs("bpairs"); 1976 node_pairs.insert("1,2", "1,3", "3,4", "4,1"); 1977 Model<> Mtest("Mtest"); 1978 var<> R_Wij("R_Wij", -1, 1); 1979 var<> Im_Wij("Im_Wij", -1, 1); 1980 var<> Wii("Wii", 0.8, 1.21); 1981 Mtest.add(R_Wij.in(node_pairs), Im_Wij.in(node_pairs), Wii.in(buses)); 1982 Constraint<> SOC("SOC"); 1983 SOC = pow(R_Wij, 2) + pow(Im_Wij, 2) - Wii.from(node_pairs)*Wii.to(node_pairs); 1984 Mtest.add(SOC.in(node_pairs) <= 0); 1985 Mtest.print(); 1986 auto g = Mtest.get_interaction_graph(); 1987 g.print(); 1988 CHECK(g.nodes.size()==12); 1989 CHECK(g.arcs.size()==4); 1990 Path p; 1991 p.nodes.push_back(g.get_node("Wii[1]")); 1992 p.nodes.push_back(g.get_node("Wii[2]")); 1993 p.nodes.push_back(g.get_node("Wii[3]")); 1994 p.nodes.push_back(g.get_node("Wii[1]")); 1995 CHECK(p.source_dest(g.get_node("Wii[1]"), g.get_node("Wii[1]"))); 1996 CHECK(p.length()==4); 1997 CHECK(p.cycle()); 1998 CHECK(p.to_str()=="{ Wii[1] , Wii[2] , Wii[3] , Wii[1] }"); 1999 } 2000 2001 TEST_CASE("testing SDP-BT"){ 2002 PowerNet grid; 2003 auto fname = string(prj_dir)+"/data_sets/Power/nesta_case9_bgm__nco.m"; 2004 try{ 2005 build_ACOPF(grid, ACRECT); 2006 } 2007 catch(invalid_argument& arg){ 2008 cout << "Error successfully caught: "<< endl; 2009 cout << arg.what() << endl; 2010 } 2011 grid.readgrid(fname); 2012 auto OPF=build_ACOPF(grid, ACRECT); 2013 double ub_solver_tol=1e-6, lb_solver_tol=1e-8, range_tol=1e-3, max_time = 200, opt_rel_tol=1e-2, opt_abs_tol=1e6; 2014 unsigned max_iter=1e3, nb_threads=1; 2015 SolverType ub_solver_type = ipopt, lb_solver_type = ipopt; 2016 auto nonlin_obj=true, current=true; 2017 auto SDP= build_SDPOPF(grid, current, nonlin_obj); 2018 auto res=OPF->run_obbt(SDP, max_time, max_iter, opt_rel_tol, opt_abs_tol, nb_threads=1, ub_solver_type, lb_solver_type, ub_solver_tol, lb_solver_tol, range_tol); 2019 auto lower_bound = SDP->get_obj_val(); 2020 auto lower_bound_init = get<3>(res); 2021 auto upper_bound = OPF->get_obj_val(); 2022 auto gap_init = 100*(upper_bound - lower_bound_init)/std::abs(upper_bound); 2023 auto final_gap = 100*(upper_bound - lower_bound)/std::abs(upper_bound); 2024 CHECK(gap_init>10); 2025 CHECK(final_gap<1); 2026 } 2027 2028 TEST_CASE("Second Derivative Of One Constraint Equals First Derivative Of Another") { 2029 Model<> mod; 2030 2031 var<double> t1("t1", -1.0, 1.0); mod.add_var(t1); t1.initialize_all(0.0); 2032 var<double> t2("t2", -1.0, 1.0); mod.add_var(t2); t2.initialize_all(0.0); 2033 2034 Constraint<double> c1("c1"); 2035 c1 += cos(t1 - t2); 2036 c1 == 0.0; 2037 mod.add_constraint(c1); 2038 2039 Constraint<double> c2("c2"); 2040 c2 += sin(t1 - t2); 2041 c2 == 0.0; 2042 mod.add_constraint(c2); 2043 2044 mod.min(t1 + t2); 2045 2046 mod.print(); 2047 2048 solver<> s(mod, ipopt); 2049 s.run(); 2050 2051 mod.print(); 2052 } 2053 2054 #ifdef USE_MPI 2055 TEST_CASE("testing MPI") { 2056 DebugOn("testing MPI" << endl); 2057 int worker_id; 2058 auto err_rank = MPI_Comm_rank(MPI_COMM_WORLD, &worker_id); 2059 unsigned nb_threads = 2; 2060 double tol = 1e-6; 2061 PowerNet grid1,grid2; 2062 string fname = string(prj_dir)+"/data_sets/Power/nesta_case5_pjm.m"; 2063 grid1.readgrid(fname); 2064 fname = string(prj_dir)+"/data_sets/Power/nesta_case14_ieee.m"; 2065 grid2.readgrid(fname); 2066 auto ACOPF1 = build_ACOPF(grid1,ACRECT); 2067 auto SOCOPF1 = build_SDPOPF(grid1); 2068 auto ACOPF2 = build_ACOPF(grid2,ACPOL); 2069 auto SOCOPF2 = build_SDPOPF(grid2); 2070 auto models = {ACOPF1, SOCOPF1, ACOPF2, SOCOPF2}; 2071 /* run in parallel */ 2072 run_MPI(models, ipopt, tol=1e-6, nb_threads=2); 2073 if(worker_id==0){ 2074 ACOPF1->print_solution(); 2075 ACOPF2->print_solution(); 2076 SOCOPF2->print_solution(); 2077 } 2078 MPI_Finalize(); 2079 } 2080 #endif 2081 2082