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