option display_precision 9; # model unit_commit_lin.mod # data DATA/toy.dat # data DATA/u10t24.dat # data DATA/u20t24.dat # data DATA/u30t24.dat # data DATA/u40t24.dat # data DATA/u50t24.dat # data DATA/u60t24.dat # data DATA/u70t24.dat # data DATA/u80t24.dat # data DATA/u90t24.dat # data DATA/u100t24.dat model run.mod data run.dat option solver cplexamp; option show_stats 1; #### option cplex_options 'timing 1 mipdisplay 1 threads 1 endbasis lin.bas startbasis lin.bas'; option cplex_options 'timing 1 mipdisplay 1 threads 1'; # option cplex_options 'threads=1 mipgap=1.e-12 integrality=1.e-9 absmipgap=1.e-12 qcpconvergetol=1.e-12'; # rel_boundtol=1.e-6'; # option solver gurobi_ampl; # option gurobi_options 'threads 1'; param epsilon := 1.e-6; param error; param lincost; param quadcost; param nCUTS; set TMP ordered; # ordered for being able to get index of its elements let {u in UNITS} nBRK[u] := 2; repeat { solve; #for {u in UNITS} { # for {i in 1..nBRK[u]} { # printf "%g\t", p_med[u,i]; # } # printf "\n"; #} printf "\n\nCURRENT TOTAL ABS ERROR: %.12g, CURRENT FUEL COST %.12g\n", sum {u in UNITS, t in 1..T} abs(a[u] * y[u,t] + b[u] * p[u,t] + c[u] * p[u,t] * p[u,t] - fuel[u,t]), sum {u in UNITS, t in 1..T} (a[u] * y[u,t] + b[u] * p[u,t] + c[u] * p[u,t] * p[u,t]); let nCUTS := 0; for {u in UNITS} { printf "\n\nUNIT %s, CURRENT ABS ERROR %.12g, CURRENT FUEL COST %.12g\n", u, sum {t in 1..T} abs(a[u] * y[u,t] + b[u] * p[u,t] + c[u] * p[u,t] * p[u,t] - fuel[u,t]), sum {t in 1..T} (a[u] * y[u,t] + b[u] * p[u,t] + c[u] * p[u,t] * p[u,t]); let TMP := setof {i in 1..nBRK[u]} p_med[u,i]; printf "initial TMP:"; display TMP; for {t in 1..T : p[u,t] > 0 and p[u,t] not in TMP} { let quadcost := a[u] * y[u,t] + b[u] * p[u,t] + c[u] * p[u,t] * p[u,t]; let lincost := max {i in 1..nBRK[u]} ((alpha_med[u,i] - beta_med[u,i] * p_med[u,i]) * y[u,t] + beta_med[u,i] * p[u,t]); ### let error := abs(quadcost - fuel[u,t])/quadcost; let error := abs(quadcost - lincost)/quadcost; if error > epsilon then { printf "inserting %.15g\n", p[u,t]; printf "quadratic costs: %.12g\n", quadcost; printf "linearized costs: %.12g (%.12g))\n", fuel[u,t], lincost; printf "error: %.12g\n", error; let TMP := TMP union {p[u,t]}; let nCUTS := nCUTS + 1; let nBRK[u] := card(TMP); # display {i in 1..nBRK[u]} member(i,TMP); let {i in 1..nBRK[u]} p_med[u,i] := member(i,TMP); # display {i in 1..nBRK[u]} p_med[u,i]; } } printf "current linear cost constraints at:\n"; display {i in 1..nBRK[u]} p_med[u,i]; } display sum{u in UNITS} nBRK[u], nCUTS, cost; } while nCUTS > 0; display p; display on; display off; display y; display cost;