Originally posted by: ALessin

I am having difficulty formulating the following constraint:

where , I but can't figure out the best way to implement it. 

I have been trying: 

forall (i in Nodes)

      ctSAMCoverage:
          (sum(t in Type, j in Nodes: (j!=i)) SAMProb[i][j][t]*Locate[j][t])*(sum(t in Type) Locate[i][t]) >= sum(t in Type) SAMCov[t]*Locate[i][t];     

but get errors saying "Implicit slicing not possible for "n" and "Circular dependency declaring "SAMProb"." Any help would be greatly appreciated. 

Thank you,

Aaron

Hi,

can you post your .mod and .dat ?

I wrote

 using CP;
 
 range Nodes=1..4;
 range Type=1..2;
 
 dvar int Locate[Nodes][Type];
 int SAMCov[i in Type]=i;
 float SAMProb[i in Nodes][j in Nodes][t in Type]=i*j*t;
 
 
 subject to
 {
 forall (i in Nodes)

      ctSAMCoverage:
          (sum(t in Type, j in Nodes:
          (j!=i)) SAMProb[i][j][t]*Locate[j][t])*(sum(t in Type) Locate[i][t]) >= sum(t in Type)
          SAMCov[t]*Locate[i][t];    
          
          
        
      }          

and that worked fine

regards

Originally posted by: ALessin

Alex,

Attached are my .mod and .dat files. I now get a solution to the original formulation, however, when I fix the solution values and try to solve for the dual constraint values, I get the error: "No dual solution exists" on line 129. Perhaps you will be able to tell why or suggest an alternative way to code the above constraint in question. Thank you for your help.

Hi,

the way you wrote ctSAMCoverage: leads to QP

So what I would do is a variable change:

dvar float changeVar[n in Nodes][m in Nodes][t in Type];

and then

forall(n,m in Nodes,t in Type) (Locate[m][t]==0) => (changeVar[n][m][t]==0);

forall(n,m in Nodes,t in Type) (Locate[m][t]==1) => (changeVar[n][m][t]==SAMProb[n][m][t]*sum(t in Type) Locate[n][t]);

and then I would change ctSAMCoverage: into

forall (n in Nodes)
      ctSAMCoverage:
          (sum(t in Type, m in Nodes: (m!=n)) changeVar[n,m,t] )>= sum(t in Type) SAMCov[t]*Locate[n][t];          
       

and then your model works fine.

I tweeted https://twitter.com/AlexFleischer1/status/843806521918132225

for this exact use case

regards

NB:

Your complete .mod

/*********************************************
 * OPL 12.7.0.0 Model - 3x3 Example Code
 * Author: Lessin
 * Creation Date: Feb 22, 2017 at 12:52:27 PM
 *********************************************/
int NumSAMTypes = ...;   // Number of SAM Types
range Type = 1..NumSAMTypes;
int NumSAMNodes = ...;   // Number of nodes
range Nodes = 1..NumSAMNodes;
int NumArcNodes = ...;   // Number of arcs
range DualVar = 0..NumArcNodes+1;
int NumPrefPts = ...;    // Number of preferential coverage points
range Prefpts = 1..NumPrefPts;

 

int MaxLocate[Type] = ...;
float ExpTypeWeight[Type] = ...;
float PrefCov[Prefpts] = ...;
float SAMCov[Type] = ...;

// Create a record to hold information about each arc
tuple arc {
   int fromnode;
   int tonode;
}

// Get the set of arcs

{arc} Arcs = ...;

float EdgeWeights1[Nodes][Arcs] = ...;
float EdgeWeights2[Nodes][Arcs] = ...;
float EdgeWeights3[Nodes][Arcs] = ...;
float EdgeWeights[n in Nodes][a in Arcs][t in Type]=(t==1)?(EdgeWeights1[n][a]):((t==2)?(EdgeWeights2[n][a]):(EdgeWeights3[n][a]));

float PrefProb1[Prefpts][Nodes] = ...;
float PrefProb2[Prefpts][Nodes] = ...;
float PrefProb3[Prefpts][Nodes] = ...;
float PrefProb[p in Prefpts][n in Nodes][t in Type]=(t==1)?(PrefProb1[p][n]):((t==2)?(PrefProb2[p][n]):(PrefProb3[p][n]));

float SAMProb1[Nodes][Nodes] = ...;
float SAMProb2[Nodes][Nodes] = ...;
float SAMProb3[Nodes][Nodes] = ...;
float SAMProb[n in Nodes][m in Nodes][t in Type]=(t==1)?(SAMProb1[n][n]):((t==2)?(SAMProb2[n][n]):(SAMProb3[n][n]));

// The network flow model has decision variables indexed on
// the arcs.

dvar float+ z;
dvar boolean Locate[Nodes][Type];
dvar float Pi[d in DualVar];
//float Duals[Arcs];
int Locate3[i in Nodes][t in Type];

dvar float changeVar[n in Nodes][m in Nodes][t in Type];

maximize z;
subject to {

 

forall(n,m in Nodes,t in Type) (Locate[m][t]==0) => (changeVar[n][m][t]==0);

forall(n,m in Nodes,t in Type) (Locate[m][t]==1) => (changeVar[n][m][t]==SAMProb[n][m][t]*sum(t in Type) Locate[n][t]);

 

    //Limit number of located sites for each SAM battery type
    forall (t in Type)
      ctNodeLocation:
          sum(n in Nodes) Locate[n][t] <= MaxLocate[t];
          
    //Only allow one SAM battery at a given location
    forall (n in Nodes)
      ctLocationLimitation:
          sum(t in Type) Locate[n][t] <= 1;
          
    //Ensure coverage of preferential coverage points
    forall (p in Prefpts)
      ctPrefPointCoverage:
          sum(t in Type, n in Nodes) PrefProb[p][n][t]*Locate[n][t] >= PrefCov[p];
          
    //Ensure coverage of SAM batteries
//    forall (n in Nodes)
//      ctSAMCoverage:
//          (sum(t in Type, m in Nodes: (m!=n)) SAMProb[n][m][t]*Locate[m][t])*(sum(t in Type) Locate[n][t]) >= sum(t in Type) SAMCov[t]*Locate[n][t];          
    
    
    forall (n in Nodes)
      ctSAMCoverage:
          (sum(t in Type, m in Nodes: (m!=n)) changeVar[n,m,t] )>= sum(t in Type) SAMCov[t]*Locate[n][t];          
        
    
          
    Pi[NumArcNodes+1]-Pi[0] >= z;
    Pi[0] == 0;
 
   // Dual constraints
     forall (a in Arcs)
     ctDual:
     Pi[a.tonode]-Pi[a.fromnode] <= sum (t in Type, n in Nodes) ExpTypeWeight[t]*EdgeWeights[n][a][t]*Locate[n][t];
}

execute DISPLAY {
writeln("\nLocate SAM batteries at:\n");
   for(var t in Type)
   for (var n in Nodes)
      if(Locate[n][t] > 0)
        writeln("SAM type ",t," at Node ",n);
   }
   
 
main
{
var status = 0;
thisOplModel.generate();
//cplex.epgap=0.015; //Ability to set the optimality gap to speed up solution time
if (cplex.solve())
{ //Solves the original problem
thisOplModel.postProcess();
writeln("\nInteger Model");
writeln("OBJECTIVE: ",cplex.getObjValue()); //Displays objective value of original problem
for(var a in thisOplModel.Arcs)writeln("dual CT value for:",a,"= ",thisOplModel.ctDual[a].dual); //Displays the "undefined" dual constraint values
write("\nSAMSolution = [") //Displays the SAM battery solution vector to be used for plotting purposes in MATLAB
for(var t in thisOplModel.Type) for (var n in thisOplModel.Nodes) if(thisOplModel.Locate[n][t] > 0) write(n," ")
write("];\n");
}

//Fixes the values of the original solution (Locate) in the relaxed problem (Locate3)
var Locate2=new Array(thisOplModel.Nodes.size);
for (var n in thisOplModel.Nodes) Locate2[n]=new Array(thisOplModel.Type.size);
for(var n in thisOplModel.Nodes) for(var t in thisOplModel.Type) Locate2[n][t]=thisOplModel.Locate[n][t];
for(var n in thisOplModel.Nodes) for(var t in thisOplModel.Type) thisOplModel.Locate3[n][t]=Locate2[n][t];

writeln("\nfrozen =",thisOplModel.Locate3); //Displays the fixed solution array from the original problem

//Sets the upper and lower bounds of the relaxed solution to the solution of the original problem
for(var n in thisOplModel.Nodes) for(var t in thisOplModel.Type) thisOplModel.Locate[n][t].UB=thisOplModel.Locate3[n][t];
for(var n in thisOplModel.Nodes) for(var t in thisOplModel.Type) thisOplModel.Locate[n][t].LB=thisOplModel.Locate3[n][t];
thisOplModel.convertAllIntVars();
if (cplex.solve())
{//Solves the relaxed problem
thisOplModel.postProcess();
writeln("\nRelaxed Model");
writeln("OBJECTIVE: ",cplex.getObjValue()); //Displays the objective value of the relaxed problem (should match original)
for(var a in thisOplModel.Arcs)writeln("dual CT value for:",a,"= ",thisOplModel.ctDual[a].dual); //Displays the dual values which represent the optimal path the enemy use to traverse the region
write("\nPathSolution = [ 0 ") //Displays the path solution vector to be used for plotting purposes in MATLAB
for(var a in thisOplModel.Arcs) if(thisOplModel.ctDual[a].dual > 0) write(a.tonode," ");
write("];")
write("\nTypesSolution = [") //Displays the types of SAM battery solution vector to be used for plotting purposes in MATLAB
for(var t in thisOplModel.Type) for (var n in thisOplModel.Nodes) if(thisOplModel.Locate[n][t] > 0) write(t," ");
write("];")
}
}
 

 

Originally posted by: ALessin

Thank you very much, Alex! The model is working now, except I'm getting "undefined" values for the ctDual values I solve for at the end of my code. Are the above changes somehow not incorporated correctly when I fix the original solution values and solve the relaxed problem? Thank you for your continued help!

Hi,

at https://www.ibm.com/developerworks/community/forums/html/threadTopic?id=aa9aa3db-4fbc-4209-a767-5b5e54902cbd&ps=25

I gave many options in order to write x*b

Let's use the linearize way.

Let us turn

forall(n,m in Nodes,t in Type) (Locate[m][t]==0) => (changeVar[n][m][t]==0);

forall(n,m in Nodes,t in Type) (Locate[m][t]==1) => (changeVar[n][m][t]==SAMProb[n][m][t]*sum(t in Type) Locate[n][t]);

into

forall(n,m in Nodes,t in Type)
  {
 
   0*Locate[m][t]<=changeVar[n][m][t];
changeVar[n][m][t]<=card(Type)*SAMProb[n][m][t]*Locate[m][t];

changeVar[n][m][t]<=SAMProb[n][m][t]*sum(t in Type) Locate[n][t]-0*(1-Locate[m][t]);
changeVar[n][m][t]>=SAMProb[n][m][t]*sum(t in Type) Locate[n][t]-card(Type)*SAMProb[n][m][t]*(1-Locate[m][t]);  
  }

then in the results you will see for example

dual CT value for: <25 16>= 0
dual CT value for: <16 25>= 0
dual CT value for: <26 17>= 0
dual CT value for: <17 26>= 0
dual CT value for: <27 18>= 0
dual CT value for: <18 27>= 1
dual CT value for: <19 29>= 0
dual CT value for: <29 19>= 0
dual CT value for: <20 30>= 0
dual CT value for: <30 20>= 0
dual CT value for: <21 31>= 0

regards

Originally posted by: ALessin

Thank you very much for all of your help, Alex! It works perfectly now.