Originally posted by: felycite28

Hi folks,

I have

parameter: a,b,x

variable c

How can I linearize:

a*e^{-bc}+b*e^{-ac}<=x  for CPLEX ? Is there any special function in CPLEX?

Thank you in advance

Hi,

why not using a piecewise linear function as an approximation ?

See

https://www.ibm.com/developerworks/community/forums/html/topic?id=dd4fa84b-4bcb-4be4-8887-0b55e8c8f7f2&ps=25

You could also try to use CP

regards

Originally posted by: felycite28

Hello Alex ,

Thank you for the answer but I could not understand it well since the function is not shown. Is there any other example which I can use for my exponential function above? "ln (a)" or does not work on CPLEX right ?

Even it works, it would not solve my function anyway but if you have any other ideas I will appreciate  ,I really stuck here .

Thank you in advance

Hi,

with the link I gave you you may approximate any function with a piecewise linear function.

Is that ok with you ? Or you are ok to use CPO ?

regards

Originally posted by: felycite28

Hello Alex,

It will be ok, I will apply it ,  but in the link you sent I could not understand the function which we need to make linear approximation , I mean original function ?

Can you tell me the function itself ?

Thank you in advance

Hi,

let me give you an example.

Suppose you need to linearize x*x+log(x)

int sampleSize=10000;
float s=1;
float e=10;

float x[i in 0..sampleSize]=s+(e-s)*i/sampleSize;

int nbSegments=20;

float x2[i in 0..nbSegments]=(s)+(e-s)*i/nbSegments;
float y2[i in 0..nbSegments]=x2[i]*x2[i]+log(x2[i]);  // Here is the function

float firstSlope=0;
 float lastSlope=0;
 
 tuple breakpoint // y=f(x)
 {
  key float x;
  float y;
 }
 
 sorted { breakpoint } breakpoints={<x2[i],y2[i]> | i in 0..nbSegments};
 
 float slopesBeforeBreakpoint[b in breakpoints]=
 (b.x==first(breakpoints).x)
 ?firstSlope
 :(b.y-prev(breakpoints,b).y)/(b.x-prev(breakpoints,b).x);
 
 pwlFunction f=piecewise(b in breakpoints)
 { slopesBeforeBreakpoint[b]->b.x; lastSlope } (first(breakpoints).x, first(breakpoints).y);
 
 assert forall(b in breakpoints) f(b.x)==b.y;
 
 float maxError=max (i in 0..sampleSize) abs(x[i]*x[i+log(x[i])-f(x[i]));
 float averageError=1/(sampleSize+1)*sum (i in 0..sampleSize) abs(x[i]*x[i+log(x[i])-f(x[i]));
 
 execute
{
writeln("maxError = ",maxError);
writeln("averageError = ",averageError);
}

will do the job and f will be the piecewise linear approximation

regards

https://www.linkedin.com/pulse/what-optimization-how-can-help-you-do-more-less-zoo-buses-fleischer/

Originally posted by: felycite28

Hello Alex,

Million thanks again . You taught me  a lot during my studies.  Now , I applied pwl function for one of my model . In the second model I have a constraint like that :

f(x)=1-e^(a-b) here a and b are decision variables .

1) Can I apply pwl for  it again ?

2) If yes , how can I express e (euler) and power in cplex as you expressed in previous example ?

3) For execution part , which I did not yet , how will i insert this approximation part to the model ? ( This question might be silly sorry but I have never done it before , just learning)

Let me explain with exo. My math model :

Objective function

Constraint 1

Constraint 2

Constraint 3

Constraint 4 // the constraint needs to be approximated , we did it

Constraint 5

Constraint 6 so on
 

Shall I directly copy and paste this approximation job into the original model ?

Million thanks

Hi,

you do the same but with your function:

 

    int sampleSize=10000;
    float s=0;
    float e=5;

    float x[i in 0..sampleSize]=s+(e-s)*i/sampleSize;

    int nbSegments=100;

    float x2[i in 0..nbSegments]=(s)+(e-s)*i/nbSegments;
    float y2[i in 0..nbSegments]=1-exp(x2[i]);  // Here is the function

    float firstSlope=0;
     float lastSlope=0;
     
     tuple breakpoint // y=f(x)
     {
      key float x;
      float y;
     }
     
     sorted { breakpoint } breakpoints={<x2[i],y2[i]> | i in 0..nbSegments};
     
     float slopesBeforeBreakpoint[b in breakpoints]=
     (b.x==first(breakpoints).x)
     ?firstSlope
     :(b.y-prev(breakpoints,b).y)/(b.x-prev(breakpoints,b).x);
     
     pwlFunction f=piecewise(b in breakpoints)
     { slopesBeforeBreakpoint[b]->b.x; lastSlope } (first(breakpoints).x, first(breakpoints).y);
     
     assert forall(b in breakpoints) f(b.x)==b.y;
     
     float maxError=max (i in 0..sampleSize) abs(1-exp(x[i])-f(x[i]));
     float averageError=1/(sampleSize+1)*sum (i in 0..sampleSize) abs(1-exp(x[i])-f(x[i]));
     
     execute
    {
    writeln("maxError = ",maxError);
    writeln("averageError = ",averageError);
    }
    
    
dvar float a in 0..4;
dvar float b in 0..4;

 subject to
 {
 abs(f(a-b)+10)<=0.001;
 
 }  
 
 float r=1-exp(a-b);
 execute
 {
 writeln("1-exp(a-b)=",r);
 }

regards

Originally posted by: felycite28

Hello Alex,

Thank you for the answer, I have a question for application .

Now my function is actually = f(x)=lambda *(1-e^(a[t]-b[t]) )

I had already defined a[t] and b[t]  like that

dvar float a[t]

dvar float b[t]

before objective function . Do I need to define them again lter as you did

dvar float a in 0..4;
dvar float b in 0..4;

I am really struglling to insert this approximation to my model ?  Should we open a new .mod and copy this approximation ?  or if i copy and past , I am getting many errors whe I try it. For example it says"
Variable de décision (ou expression) "a" non autorisée. 

Thank you so much

Then you could write

float lambda=1.2;
int sampleSize=10000;
    float s=0;
    float e=5;

    float x[i in 0..sampleSize]=s+(e-s)*i/sampleSize;

    int nbSegments=100;

    float x2[i in 0..nbSegments]=(s)+(e-s)*i/nbSegments;
    float y2[i in 0..nbSegments]=lambda*(1-exp(x2[i]));  // Here is the function

    float firstSlope=0;
     float lastSlope=0;
     
     tuple breakpoint // y=f(x)
     {
      key float x;
      float y;
     }
     
     sorted { breakpoint } breakpoints={<x2[i],y2[i]> | i in 0..nbSegments};
     
     float slopesBeforeBreakpoint[b in breakpoints]=
     (b.x==first(breakpoints).x)
     ?firstSlope
     :(b.y-prev(breakpoints,b).y)/(b.x-prev(breakpoints,b).x);
     
     pwlFunction f=piecewise(b in breakpoints)
     { slopesBeforeBreakpoint[b]->b.x; lastSlope } (first(breakpoints).x, first(breakpoints).y);
     
     assert forall(b in breakpoints) f(b.x)==b.y;
     
     float maxError=max (i in 0..sampleSize) abs(lambda*(1-exp(x[i]))-f(x[i]));
     float averageError=1/(sampleSize+1)*sum (i in 0..sampleSize) abs(lambda*(1-exp(x[i]))-f(x[i]));
     
     execute
    {
    writeln("maxError = ",maxError);
    writeln("averageError = ",averageError);
    }
    
range t=1..5;    
dvar float a[t] in 0..4;
dvar float b[t] in 0..4;

 subject to
 {
 forall(i in t) abs(f(a[i]-b[i])+10)<=0.001;
 forall(i in t:i+1 in t) a[i]<=a[i+1]-0.1;
 }  
 
 float r[i in t]=1-exp(a[i]-b[i]);
 execute
 {
 writeln("1-exp(a-b)=",r);
 }

regards

Originally posted by: felycite28

Hello Alex,

I understood that part. My question is hw will I combine this approximation and original model? For exaple when we code on python , we can open a .py and code a function and open another .py it becomes main code and we call 1.py ? Like this fashion , lets say on the original.mod I have my original code, on the approxim.mod I have approximation code on cplex ; Can I call it? Or how will i integrate both?

range t=1..5;    
dvar float a[t] in 0..4;
dvar float b[t] in 0..4;  subject to
 {
 forall(i in t) abs(f(a[i]-b[i])+10)<=0.001;
 forall(i in t:i+1 in t) a[i]<=a[i+1]-0.1;
 }  

This part i your original code ?  If so , you code approxim thing before the original code ? I could not understand just combining both ?

Then you could use include

in pwlexp.mod

float lambda=1.2;
int sampleSize=10000;
    float s=0;
    float e=5;

    float x[i in 0..sampleSize]=s+(e-s)*i/sampleSize;

    int nbSegments=100;

    float x2[i in 0..nbSegments]=(s)+(e-s)*i/nbSegments;
    float y2[i in 0..nbSegments]=lambda*(1-exp(x2[i]));  // Here is the function

    float firstSlope=0;
     float lastSlope=0;
     
     tuple breakpoint // y=f(x)
     {
      key float x;
      float y;
     }
     
     sorted { breakpoint } breakpoints={<x2[i],y2[i]> | i in 0..nbSegments};
     
     float slopesBeforeBreakpoint[b in breakpoints]=
     (b.x==first(breakpoints).x)
     ?firstSlope
     :(b.y-prev(breakpoints,b).y)/(b.x-prev(breakpoints,b).x);
     
     pwlFunction f=piecewise(b in breakpoints)
     { slopesBeforeBreakpoint[b]->b.x; lastSlope } (first(breakpoints).x, first(breakpoints).y);
     
     assert forall(b in breakpoints) f(b.x)==b.y;
     
     float maxError=max (i in 0..sampleSize) abs(lambda*(1-exp(x[i]))-f(x[i]));
     float averageError=1/(sampleSize+1)*sum (i in 0..sampleSize) abs(lambda*(1-exp(x[i]))-f(x[i]));
     
     execute
    {
    writeln("maxError = ",maxError);
    writeln("averageError = ",averageError);
    }

and then you can write

include "pwlexp.mod";
    
range t=1..5;    
dvar float a[t] in 0..4;
dvar float b[t] in 0..4;

 subject to
 {
 forall(i in t) abs(f(a[i]-b[i])+10)<=0.001;
 forall(i in t:i+1 in t) a[i]<=a[i+1]-0.1;
 }  
 
 float r[i in t]=1-exp(a[i]-b[i]);
 execute
 {
 writeln("1-exp(a-b)=",r);
 }

regards

Originally posted by: felycite28

Hello Alex ,

Can I ask one more question ?

forall(i in t) abs(f(a[i]-b[i])+10)<=0.001;
 forall(i in t:i+1 in t) a[i]<=a[i+1]-0.1;

This constraints are the random constraints to  show a global.mod or are they piecewise related constrainsts? I could not understand them , sorry .

2) Actually constraint was like f(x)=lambda *(1-e^(a[t]-b[t]) ) <=5 , in this case how will I define it in the global model.  you had written

 float r[i in t]=1-exp(a[i]-b[i]);

without less than version . If constraint is as above , how should I define it ?

Thank you in advance.

Hi,

forall(i in t) abs(f(a[i]-b[i])+10)<=0.001; is a way to relax a bit the equality since the pwl is an approximation
 forall(i in t:i+1 in t) a[i]<=a[i+1]-0.1; was just to give an example of different values for a[i] you may forget that

float r[i in t]=1-exp(a[i]-b[i]);  is a way to compute the results

regards

PS:

Not many days left for https://www.linkedin.com/feed/update/urn:li:activity:6409320254585004032

Originally posted by: felycite28

Hello Alex ,

I got it I guess, we do not have to use "abs" constraint we can use any other type constraint to relax the model , which will get the model feasible solution , right ? Just I wanted to be sure if I got it correct.

Can we plot the graph of pw function on cplex ? If yes , how? Thank you so much

Hi,

you could do that through a call to python

See

https://www.ibm.com/developerworks/community/forums/html/topic?id=9361b722-1efc-43f1-9251-6e6bb4f38cf0&ps=25

or for your instance:

execute
{

// turn an OPL array into a python list
function getPythonListOfArray(_array)
{

var quote="\"";
var nextline="\\\n";

var res="[";
for(var i in _array)
{
var value=_array[i];

if (typeof(value)=="string") res+=quote;
res+=value;
if (typeof(value)=="string") res+=quote;
res+=",";
res+=nextline;
}
res+="]";
return res;
}

// Display a function with points with x and y arrays of x and y
function displayXY(x,y,pythonpath,pythonfile)
{
writeln("displayXY ",x," ",y," ",pythonpath," ",pythonfile);

var python=new IloOplOutputFile(pythonfile);
python.writeln("import matplotlib.pyplot as plt");
python.writeln("x = ",getPythonListOfArray(x))
python.writeln("y = ",getPythonListOfArray(y))
python.writeln("plt.plot(x, y)");
python.writeln("plt.xlabel('x - axis')");
python.writeln("plt.ylabel('y - axis')");
python.writeln("plt.title('xy graph')");
python.writeln("plt.show()");
python.close();
IloOplExec(pythonpath+" "+ pythonfile,true);        
}

 

 

 

}

float s=1;
float e=10;

 

int nbSegments=20;

float x2[i in 0..nbSegments]=(s)+(e-s)*i/nbSegments;
float y2[i in 0..nbSegments]=x2[i]*x2[i]+log(x2[i]);  // Here is the function

string pythonpath="C:\\Users\\IBM_ADMIN\\AppData\\Local\\Programs\\Python\\Python36\\python.exe";
string pythonfile="C:\\displayXY.py";
execute
{

displayXY(x2,y2,pythonpath,pythonfile);
}

which gives

and if you turn nbsegments to 3 you ll get

regards

Originally posted by: felycite28

Hello Alex,

I need to ask a question . can we use f function in objective function of global model? For example you have constraint  

forall(i in t) abs(f(a[i]-b[i])+10)<=0.001; you call f function from approximation.mod and it works accordingly, it is ok for me too . But technic. speaking, cam we do something like that and does it work as well ? For example,

Min z= a+b+ f(a[i]-b[i])) ; is it possible to call linearization function for objectve function as we can do for constraints ? I tried it when the number of segments are very low, it works but when the number of segments increases it can not find a solution for a long CPU time (when number of segments=5) is it normal, how can I deal with it?

Thank you so much

Hi,

well don t forget that by having many segments you increase the size of the model.

Plus you could turn a feasible model into not feasible if you do not allow enough tolerance

regards

Originally posted by: felycite28

Hello Alex,

I need to ask a question . We had discussed the pw linearization before . I tested my model with LINGO ,nonlinear solver, but as you can guess it did ot reach global optimum.

I found the link of Paul Rubin about the topic.

https://www.quora.com/How-can-I-solve-a-constrained-optimization-problem-where-the-objective-function-is-sum-of-a-linear-expression-and-an-exponential-of-linear-expression

My objective function is the same with the example of Paul.  Here it seems that we have to lienarize inside of the exp which is w. I tried to do pw linearization with the idea of SOS2 but I could not achieve . Can you show me an example how to linearize that objective function ? Another problem for me , after linearization , while solving on cplex how will I modify the objective fuction according to linearization ? Sorry, I am so new for nonlinear problems and linearization studies. Since the summer, it is still challange for me Any of the your help will save me

Thank you so much in advance

Originally posted by: felycite28

Hello Aex,

I need to ask a quick question: In the approximation model, starting and ending points are like this:

float s=1

float e=10

Can we define the interval with negative values ? Like

float s=-10

float e=10

and when I tried the negative version why it takes too much time , the model can not reach a solution after one hour ?

Million thanks in advance

Hi

 

    execute
    {

    // turn an OPL array into a python list
    function getPythonListOfArray(_array)
    {

    var quote="\"";
    var nextline="\\\n";

    var res="[";
    for(var i in _array)
    {
    var value=_array[i];

    if (typeof(value)=="string") res+=quote;
    res+=value;
    if (typeof(value)=="string") res+=quote;
    res+=",";
    res+=nextline;
    }
    res+="]";
    return res;
    }

    // Display a function with points with x and y arrays of x and y
    function displayXY(x,y,pythonpath,pythonfile)
    {
    writeln("displayXY ",x," ",y," ",pythonpath," ",pythonfile);

    var python=new IloOplOutputFile(pythonfile);
    python.writeln("import matplotlib.pyplot as plt");
    python.writeln("x = ",getPythonListOfArray(x))
    python.writeln("y = ",getPythonListOfArray(y))
    python.writeln("plt.plot(x, y)");
    python.writeln("plt.xlabel('x - axis')");
    python.writeln("plt.ylabel('y - axis')");
    python.writeln("plt.title('xy graph')");
    python.writeln("plt.show()");
    python.close();
    IloOplExec(pythonpath+" "+ pythonfile,true);        
    }

     

     

     

    }

    float s=-10;
    float e=10;

     

    int nbSegments=20;

    float x2[i in 0..nbSegments]=(s)+(e-s)*i/nbSegments;
    float y2[i in 0..nbSegments]=x2[i]*x2[i];  // Here is the function

    string pythonpath="C:\\Users\\IBM_ADMIN\\AppData\\Local\\Programs\\Python\\Python36\\python.exe";
    string pythonfile="C:\\displayXY.py";
    execute
    {

    displayXY(x2,y2,pythonpath,pythonfile);
    }

gives

But remember that you cannot compute the log of a negative value

regards

Originally posted by: felycite28

Thank you Alex,

yes yes I know for log function you are right. I had asked this question for general cases. For example my function f(x)=exp(x) where x can take negative values according to the mathematical model.

In that case I have to define

float s=-2

float e=2 the firts thing , can we define negative starting point like this? I tried this way , but it took too much time to reach the solution? Am I doing something wrong? or isn't it possible to take into account negative part of the function somehow ?

Million thanks

Hi

/*********************************************
 * OPL 12.8.0.0 Model
 * Author: AlexFleischer
 * Creation Date: Jul 12, 2018 at 11:53:32 AM
 *********************************************/

    float lambda=1.2;
    int sampleSize=10000;
        float s=-2;
        float e=2;

        float x[i in 0..sampleSize]=s+(e-s)*i/sampleSize;

        int nbSegments=100;

        float x2[i in 0..nbSegments]=(s)+(e-s)*i/nbSegments;
        float y2[i in 0..nbSegments]=exp(x2[i]);  // Here is the function

        float firstSlope=0;
         float lastSlope=0;
         
         tuple breakpoint // y=f(x)
         {
          key float x;
          float y;
         }
         
         sorted { breakpoint } breakpoints={<x2[i],y2[i]> | i in 0..nbSegments};
         
         float slopesBeforeBreakpoint[b in breakpoints]=
         (b.x==first(breakpoints).x)
         ?firstSlope
         :(b.y-prev(breakpoints,b).y)/(b.x-prev(breakpoints,b).x);
         
         pwlFunction f=piecewise(b in breakpoints)
         { slopesBeforeBreakpoint[b]->b.x; lastSlope } (first(breakpoints).x, first(breakpoints).y);
         
         assert forall(b in breakpoints) f(b.x)==b.y;
         
         float maxError=max (i in 0..sampleSize) abs(lambda*(1-exp(x[i]))-f(x[i]));
         float averageError=1/(sampleSize+1)*sum (i in 0..sampleSize) abs(exp(x[i]-f(x[i])));
         
         execute
        {
        writeln("maxError = ",maxError);
        writeln("averageError = ",averageError);
        }
        
        execute
    {

    // turn an OPL array into a python list
    function getPythonListOfArray(_array)
    {

    var quote="\"";
    var nextline="\\\n";

    var res="[";
    for(var i in _array)
    {
    var value=_array[i];

    if (typeof(value)=="string") res+=quote;
    res+=value;
    if (typeof(value)=="string") res+=quote;
    res+=",";
    res+=nextline;
    }
    res+="]";
    return res;
    }

    // Display a function with points with x and y arrays of x and y
    function displayXY(x,y,pythonpath,pythonfile)
    {
    writeln("displayXY ",x," ",y," ",pythonpath," ",pythonfile);

    var python=new IloOplOutputFile(pythonfile);
    python.writeln("import matplotlib.pyplot as plt");
    python.writeln("x = ",getPythonListOfArray(x))
    python.writeln("y = ",getPythonListOfArray(y))
    python.writeln("plt.plot(x, y)");
    python.writeln("plt.xlabel('x - axis')");
    python.writeln("plt.ylabel('y - axis')");
    python.writeln("plt.title('xy graph')");
    python.writeln("plt.show()");
    python.close();
    IloOplExec(pythonpath+" "+ pythonfile,true);        
    }

     

     

     

    }

    

    

    string pythonpath="C:\\Users\\IBM_ADMIN\\AppData\\Local\\Programs\\Python\\Python36\\python.exe";
    string pythonfile="C:\\displayXY.py";
    execute
    {

    displayXY(x2,y2,pythonpath,pythonfile);
    }
 

gives

regards

Originally posted by: felycite28

Thank you so much

Maybe non linear my function is too complicated and it takes too much time to compute

Thank you so much again