Originally posted by: wzyryan

I think feasOpt returns a minimization result while I was using IloMaximize. The results returned by calling IloMaximize and IloMinimization are exactly the same. Did I do something wrong while calling the feasOpt functions?

IloEnv env;
    
        IloInt w=128;
        IloInt h=128;
    typedef IloArray<IloIntArray>    IntMatrix;
        IntMatrix mat(env, w);
        for(IloInt i=0;i<w;i++){
                mat[i]=IloIntArray(env,h);
        }

        ifstream inf("out.txt");
        ofstream onf("mat.txt");
        if(!inf) cout<<"error"<<endl;

    for(IloInt j=0;j<h;j++){
        for(IloInt i=0;i<w;i++){
            inf>>mat[i][j];            //图像灰度化，储存像素点的灰度信息
        }
    }
        inf.close();

        for(IloInt j=0;j<h;j++){
        for(IloInt i=0;i<w;i++){
                        if(mat[i][j] >= 127) mat[i][j] = 255;     //图像二值化，作为Fpq输入
                        else mat[i][j] = 0;

                        onf<<mat[i][j]<<" ";
        }
                onf<<endl;
    }

        typedef IloArray<IloIntVarArray> IntVarMatrix;  
        IloInt i,j;
        IloInt n=4;
        IloInt width=32;
        IloInt height=32;
        IntVarMatrix hmat(env, width); //前面定义过的变量矩阵，hmat是其的一个对象，范围是width
        for(i=0;i<width;i++){
                hmat[i]=IloIntVarArray(env,height,0,255);
        }

        typedef IloArray<IloNumExprArray> NumExprMatrix;
        NumExprMatrix fmat(env, w);
        for(i=0;i<w;i++){
                fmat[i]=IloNumExprArray(env,h);
        }
        for(i=0;i<w;i++){
                for(j=0;j<h;j++){
                        fmat[i][j]=IloNumExpr(env, 0);
                }
        }

        IloModel model(env);              //to construct a modeling object named model, within the environment env
        IloIntVarArray t(env);
        IloRangeArray c(env);

        t.add(IloIntVar(env,20,240));
        t.add(IloIntVar(env,20,240));
        c.add( t[0]-t[1]>=0 );
        model.add(c);
        for(j=0;j<height;j++){
            for(i=0;i<width;i++){    //big pixel
                //IloInt count=0;
                for(IloInt q=0;q<n;q++){
                    for(IloInt p=0;p<n;p++){    //small pixel
                       //for computing the light intensity

                       IloInt x,y,startx,starty,boundryx,boundryy;  //实现图像边缘的3个像素范围的遍历
                       if(j>=3){
                           starty=j-3;
                       }else{
                           starty=0;
                       }
                       if(i>=3){
                           startx=i-3;
                       }else{
                           startx=0;
                       }
                       if(j+3>height){
                           boundryy=height;
                       }else{
                           boundryy=j+3;
                       }
                       if(i+3>width){
                           boundryx=width;
                       }else{
                           boundryx=i+3;
                       }

                       for(y=starty;y<boundryy;y++){
                           for(x=startx;x<boundryx;x++){     //convolution range
                               if(((n*x+(float)n/2-n*i-p)*(n*x+(float)n/2-i*n-p)+(n*y+(float)n/2-j*n-q)*(n*y+(float)n/2-j*n-q))<=9*n*n){
                               
                                                                   fmat[n*i+p][n*j+q] += (hmat[x][y])*(exp(-((x+1/2-i-(float)p/n)*(x+1/2-i-(float)p/n)+(y+1/2-j-(float)q/n)*(y+1/2-j-(float)q/n))/2))/6.28;     //小像素点的光强的叠加计算
                               }
                           }
                       }
                                           
                    }
                }
               
            }
        }

             
        for(int p=0;p<w;p++){    //logic constraints
   for(int q=0;q<h;q++){

           if(mat[p][q] == 255) model.add(t[0] <= fmat[p][q]);
           else model.add(fmat[p][q] <= t[1]);
           
   }
}
        
        
        
        model.add(IloMaximize(env, t[0]-t[1]));

        IloCplex cplex(model);
        
        cplex.exportModel("model.lp");
        cplex.solve();

        IloNumArray vals(env);
    IloNumArray infeas(env);

        if ( cplex.getStatus() == IloAlgorithm::Infeasible ||
           cplex.getStatus() == IloAlgorithm::InfeasibleOrUnbounded ) {
           env.out() << endl << "*** Model is infeasible ***" << endl << endl;
        }

        // begin feasOpt analysis
        
        cplex.setOut(env.getNullStream());
        IloNumArray lb(env);
        IloNumArray ub(env);
        
        // first feasOpt call
        
        env.out() << endl << "*** First feasOpt call ***" << endl;
        env.out() << "*** Consider all constraints ***" << endl;
        int rows = c.getSize();
        lb.add(rows, 1.0);
        ub.add(rows, 1.0);
        
        if ( cplex.feasOpt(c, lb, ub) ) {
        env.out() << endl;
        cplex.getInfeasibilities(infeas,c);
        env.out() << "*** Suggested bound changes = " << infeas << endl;
        env.out() << "*** Feasible objective value would be = "
               << cplex.getObjValue() << endl;
        env.out() << "Solution status    = " << cplex.getStatus() << endl;
        env.out() << "Solution obj value = " << cplex.getObjValue() << endl;
        cplex.getValues(vals, t);
        env.out() << "Values             = " << vals << endl;
        env.out() << endl;
        }
        else {
        env.out() << "*** Could not repair the infeasibility" << endl;
        throw (-1);
        }

                // second feasOpt call
        
        env.out() << endl << "*** Second feasOpt call ***" << endl;
        env.out() << "*** Consider all but first constraint ***" << endl;
        
        lb[0]=ub[0]=0.0;
        
        if ( cplex.feasOpt(c, lb, ub) ) {
        env.out() << endl;
        cplex.getInfeasibilities(infeas,c);
        env.out() << "*** Suggested bound changes = " << infeas << endl;
        env.out() << "*** Feasible objective value would be = "
            << cplex.getObjValue() << endl;
        env.out() << "Solution status    = " << cplex.getStatus() << endl;
        env.out() << "Solution obj value = " << cplex.getObjValue() << endl;
        cplex.getValues(vals, t);
        env.out() << "Values             = " << vals << endl; 
        env.out() << endl;
        }
        else {
            env.out() << "*** Could not repair the infeasibility" << endl;
            throw (-1);
        }

        env.end();
        getchar();

feasOpt neither minimizes nor maximizes your objective function. It finds a minimal relaxation of the constraints that makes the model feasible. Whether the model is feasible or not, and what to do to make it feasible if it is infeasible as originally specified, is independent of the objective function.

Originally posted by: wzyryan

Hi Paul,

Thank you for your reply! I've exported the model. The model exported is obviously wrong since its obj is wrong. As is indicated in line 117 of my code, the obj should be t[0]-t[1]. However, the obj in the exported file is X1 - X2 + X1027. I wonder if I did something wrong while setting the obj in my code.

Regards,

Ryan

I don't think your model is necessarily wrong. If you don't explicitly set names for your variables then CPLEX will automatically assign names x1, ..., xN to your variables. Also the Concert layer always adds a fixed variable to model constant terms in the objective function, even if that constant term is 0. In your case this added variable is x1027. As you can see in the "Bounds" section of your model, this variable is fixed to 0.

Originally posted by: wzyryan

Thanks for your reply! Now, I see. The X1 and X2 in the file represents t[0] and t[1] while X3 to X1026 represents the elements of hmat[32][32], right?

regards

Ryan

Yes, that is correct. However, you should not rely on that naming scheme. If you want to look at the model in the LP file you should explicitly name all the IloNumVar instances in your program using IloNumVar::setName() or a constructor that takes a name argument. For example, if you replace

t.add(IloIntVar(env,20,240));
t.add(IloIntVar(env,20,240));

by

t.add(IloIntVar(env,20,240, "t0"));
t.add(IloIntVar(env,20,240, "t1"));

then you should see the names t0 and t1 in the LP file (I did not use 't[0]' and 't[1]' since '[' and ']' are not valid name characters in LP files and will be replaced).

Originally posted by: wzyryan

Hi Daniel，

I replaced my code with yours. Now, t0 and t1 can be seen explicitly in the LP file. However, the X3 X4 X5...X1026 changes into X3 X5 X7...X2047. I wonder how can I name those variables from X1 to X1024.

Regards

Ryan

Like I said: you have to invoke IloNumVar::setName() on the IloNumVar instances in your code.

Originally posted by: wzyryan

As you can see, I can name t0 and t1 with the following code.

t.add(IloIntVar(env,20,240, "t0"));
t.add(IloIntVar(env,20,240, "t1"));

However, there are 32*32=1024 variables in hmat[][]. It's silly to just list from x0 to x1023.  How can I add those digit suffixes correctly?

How about using a loop in conjunction with either sprintf() or a std::stringstream?

Originally posted by: wzyryan

Hi Daniel,

I've got a set of output results and I have a question about them. As you can see, the ranges of t0 and t1 are set between 20 and 240 in the code above.

When the range is set differently, corresponding results are listed below.

(0,255)          result: t0=0,   t1=0

(10,255)        result: t0=10, t1=14

(20,255)        result: t0=20, t1=28

(30,255)        result: t0=30, t1=42

(40,255)        result: t0=40, t1=55

It is weired that t0 always equals the lower bound. Is this a unique feature of feasOpt when it relax constraints?

Regards

Ryan

Remember what Paul said: feasopt uses a completely different objective.

I think t0 always ends up at its lower bound because that is the most feasible value. So to reduce constraint violation it makes sense to put it at that value.