'inf' as a constraint

[lily ye]

It's related to fuzzy thoery. If my objective function is min x+y, but one of my  constraint is : inf (a+b)<0 or inf a+b ,how to deal with this constraint?

[Johan Löfberg]

What on earth does inf (a+b)<0 even mean?

[lily ye]

it probably means its minimize vaule is less than o

[Johan Löfberg]

Probably? So you don't even know what you are trying to say. How can you expect anyone else to undserstand what you are doing then?

Are a and b vectors? If so min(a+b)<=0 (strict inequalities are not supported). If a and b are scalars, it makes no sense to use the min operator

[lily ye]

Thanks for your reply.They are both vectors,and i just saw it in a paper .Thank you very much.

[Johan Löfberg]

No. You sent me a picture with the model, and in that picture inf meant infimizer of a scalar objective in a bilevel setup. Definitely not a min operator of a vector.

[lily ye]

So...is there any difference?how to deal with it?

[Johan Löfberg]

There is a world of difference between saying min(x) <= 0 (which is a nonconvex constraints, but rather easily written using a MILP formulation), and reasoning over bilevel formulations. It's like comparing apples and moon landers...

However, the model in the picture you sent is could possibly be silly. Adding the infimizer operator inf f(x) is redundant. No need to have a constraint inf f(x) <= 0 on the optimal objective in the inner program. Simply replace it with  the feasibility condition f(x) <= 0, and if this is feasible, the clearly the infimizer is also non non-positive. I looked up the paper you are reading, and that is precisely what they are doing, and you end up with a simple SOCP

[lily ye]

they said they drop the t for simplicity，so it is a constraint (min(a+b)>0)? but it said Warning: Solver not found (cplex) 

[Johan Löfberg]

The fact that they simplify *notation* by dropping t has absolutely nothing to do with this. It’s completely uninteresting. The important thing is that the infimizer in inf f <= simply is dropped, and a obviously sufficient feasibility condition f <=0 is used instead

Well, YALMIP is pretty clearly telling you what is wrong through the warning…

[lily ye]

"which is a nonconvex constraints, but rather easily written using a MILP formulation"?what does it mean?

[Johan Löfberg]

Don't care about that, because that discussion was before we understood what you talked about

inf in your paper means infimizer in an inner program, which is removed and replaced with simple feasibility instead

[lily ye]

I find that :You have unbounded variables in an implication leading to a lousy big-M

relaxation.It has nothing to do with this constraints? it is about "implies" in my model? 

Besides, you mean i don't need to write down this constraints?caz it has been replaced by the following constraints?

[Johan Löfberg]

Why do you have implies? There are no implies in the model in the paper.

Read the paper, and remove the word inf in (35). That's your model. Well, beyond that, I think there are some quadratic constraints that you have to write manually as SOCP cones also, YALMIP will probably not do that automatiaclly.

[lily ye]

i  try to use one variables to define the charge and discharge together，so there is  "implies" to judge where it is below or above 0, and use different formulation...maybe i misunderstand the"implies"? I thought it equals to "if"?  

it warings in my "implies " line  :"Strict inequalities are not supported. A non-strict has been added instead "

I  will try to transfer the  quadratic constraints manually ,and quadratic constraint means  "variables^2"?

[lily ye]

And i also use implies to do linear piecewise to replace one quadric formulaiton 

[lily ye]

Something like implies(p(x)<0,d(x)=p(x)*1), implies(0<p(x)<1,d(x)=p(x)*2), implies(1<p(x)<2,d(x)=p(x)*3)

[Johan Löfberg]

so you are not implementing the paper you are discussing?

The complicating quadratic constraint is 35b which is a nonconvex quadratic constraint due to the z*s term, but since they z and s are non-negative, it is convex on the feasible domain and can be written as as an SOCP constrain as x^2 <= 4*z*s, z>=0, s>=0 is equivalent to norm([z-s;x]) <= z+s

the warning about strict inequalities is simply that you are using strict inequalities. That is not allowed since it doesn't make sense in continuous numerical optimization, YALMIP will interpret > as >=0 and < as <=

[Johan Löfberg]

which means you will get an absolutely horrible nonlinear nonconvex integer program...

[lily ye]

Thanks for your patience.Thanks for your help.

[lily ye]

I read the reference paper.

So you mean

sdpvar z,x,s

4zs>=x^2+constant   equals to  z+s>=sqrt(x^2+(z-s)^2+constant) ?

[Johan Löfberg]

No, I already told you to use a norm operator so you get an SOCP-representable expression

[lily ye]

how about my constant value?  can I write it as  norm([z-s;x])+constant <= z+s   ？

[lily ye]

If it is socp formulation，what kind of solver can i use？i don't know why it wanings sovler not found. I search the internet only to find that 

Mixed Integer Quadratic programming (commercial)

CPLEX (free for academia), GUROBI (free for academia), MOSEK (free for academia), XPRESS (free for academia)

I try to download the cplex,but it said

Access Denied

You don't have permission to access ....on this server.

[Johan Löfberg]

x^2+a^2 <= 4*z*s would consequently be norm(z-s;x;a]) <= z+s

[Johan Löfberg]

where a^2 corresponds to constant of course, so you get norm([z-s;x;sqrt(constant)]) <= z+s

[Johan Löfberg]

The fact that you don't manage to download cplex is an issue you will have to discuss with IBM/ILOG.

Much easier to get a license and download for mosek or gurobi though

[lily ye]

I'm sorry,I didn't notice that ,Cost is my objective! So..in this case , there is a standard SOCP in my constraints,there is no need to use a linear approximation?

Maybe I can try gurobi?  Why Warning: Solver not applicable (gurobi)  I failed...

[Johan Löfberg]

if you have socp constraints, you obviously don't have to approximate them as linear constraints if you have an socp solver

the message tells you that you don't have an socp though...

[lily ye]

is norm([z-s;x;sqrt(constant)]) <= z+s  a socp constraint ?

You mean I'm not allowed to use gurobi as a solver to solve SOCP ?

What shuold I do? to purchase it?

[Johan Löfberg]

Yes, that is clearly a very simple SOCP constraint (as long as x,z,s are linear expressions, and constant really is a constant), hence any SOCP solver is applicable

[Johan Löfberg]

what I am telling you is that you are loking at the wrong constraint or expression. That constraint is trivially an SOCP, but something else in your model is not SOCP representable

[lily ye]

I see, I need to  promise every constraint is SOCP! Professer,could you give me some advices on SOCP?I just know little about it..

[Johan Löfberg]

Well, first I would assume you are reading books on convex conic optimization as you appear to work on this, such as https://web.stanford.edu/~boyd/cvxbook/

Why don't you post (as simple as possible version) of your non-working code

[Johan Löfberg]

code does not work

Undefined function or variable 'time'.

 

[Johan Löfberg]

Code does not run

Error using sdpvar/interp1 (line 50)

Data should be vectors for scalar interpolant

 

However, your set of constraints are applicable to Gurobi so I don't know why you are talking about stuff not being possible to solve in gurobi. They are infeasible though, as G6 is infeasible

optimize(G6)

Solved in 0 iterations and 0.03 seconds

Infeasible model

ans = 

  struct with fields:

    yalmiptime: 0.2801

    solvertime: 0.0389

          info: 'Infeasible problem (GUROBI-GUROBI)'

       problem: 1

[Johan Löfberg]

also, yalmip is screaming about you to get rid of the strict constraints. they are not supported

[lily ye]

Ahh，what's wrong with that one？ 

AND I really don't understand why It warings:Warning: Solver not applicable (gurobi) 

Is it raleted to the different version? MATLAB with the modeling tool YALMIP [47] and GUROBI7.5.2 solver on a labtop with an Intel Core i7-5600 CPU 2.60 GHz and 8 GB of RAM.

[Johan Löfberg]

and this

z(t)+s(t)>norm([z(t)-s(t);tau(t);sqrt(1)])+3*variance_max(t)

is not consistent with the paper. In the paper it says 4*z*s >= tau^2 + 3*h^2*sigma. 

You have forced h to be 1 hence 4*z*s >= tau^2 + 3*sigma^2. , i.e. z + s >= norm([z-s;tau;sqrt(3)*sigma]). You have written z + s >= norm([z-s;tau]) + 3*sigma^2, which definitely isn't the same (square and expand both sides and simplify). You must of course confirm that any reformulation is correct, not just use them randomly

[lily ye]

I write help gurobi,and it does work

 gurobi  Solve an LP, QP, QCP, SOCP, or MIP using the Gurobi Optimizer

    result = gurobi(MODEL, PARAMS)

 

    The Gurobi MATLAB interface can be used to solve optimization problems

    of the form:

 

    minimize           c'*x + x'*Q*x + alpha

       x

    subject to         A*x = b,

                       l <= x <= u,

                       some xj integral,

                       some xk must lie within second order cones,

                       x'*Qc*x + q'*x <= beta,

                       some xi in SOS constraints.

 

    Many of the model components listed here are optional.  For example,

    integrality constraints may be omitted.

 

    The model structure must contain the following fields:

 

    model.A:  The linear constraint matrix. This must be a sparse matrix.

 

    model.obj: The linear objective vector (c in the above problem

               statement). You must specify one value for each column of A.

               This must be a dense vector.

 

    model.sense: The sense of the linear constraints. Allowed values are

                 '<', '=', or '>'. You must specify one value for each row

                 of A, or a single value to specify that all constraints

                 have the same sense. This must be a char array.

 

    model.rhs: The right-hand side vector for the linear constraints (b in

               the above problem statement). You must specify one value for

               each row of A. This must be a dense vector.

 

    The model structure may contain the following optional fields:

 

    model.lb: The lower bounds on the variables. When present, you must

              specify one value for each column of A. This must be a dense

              vector. When absent, each variable has a lower bound of 0.

 

    model.ub: The upper bounds on the variables. When present, you must

              specify one value for each column of A. This must be a dense

              vector. When absent, the variables have infinite upper

              bounds.

 

    model.vtype: The variable types. This char array is used to capture

                 variable integrality constraints. Allowed values are 'C'

                 (continuous), 'B' (binary), 'I' (integer), 'S'

                 (semi-continuous), or 'N' (semi-integer). When present,

                 you must specify one value for each column of A, or a

                 single value to specify that all variables should have the

                 same type. When absent, each variable is treated as

                 being continuous.

 

    model.modelsense: The optimization sense. Allowed values are 'min'

                      (minimize) or 'max' (maximize). When absent, the

                      default model sense is minimization.

 

    model.modelname: A name for the model. This is used in the log and

                     result files.

 

    model.objcon: The constant offset in the objective function (alpha in

                  the above problem statement).

 

    model.vbasis: The variable basis status vector. Used to provide an

                  advanced starting point for the simplex algorithm. You

                  would generally never concern yourself with contents of

                  this array, but would instead simply copy it from the

                  results of a previous optimization run. When present, you

                  must specify one value for each column of A. This must be

                  a dense vector.

 

    model.cbasis: The constraint status vector. Used to provide an advanced

                  starting point for the simplex algorithm. You would

                  generally never concern yourself with contents of this

                  array, but would instead simply copy it from the results

                  of a previous optimization run. When present, you must

                  specify one value for each row of A. This must be a dense

                  vector.

 

    model.Q: The quadratic objective matrix. When present, Q must be a

             square matrix whose row and column counts are equal to the

             number of columns of A. Q must be a sparse matrix.

 

    model.cones:  The second-order cone constraints. A struct array. When

                  present, each element in the array defines a single cone

                  constraint:

                      x(k)^2 >= sum(x(idx).^2),  x(k) >= 0.

                  The constraint is defined via model.cones.index = [k idx]

                  with the first entry in index corresponding to the index

                  of the variable on the left-hand side of the constraint,

                  and the remaining entries corresponding to the indices of

                  variables on the right-hand side of the constraint.

                  model.cones.index must be a dense vector.

                  For example, the two second-order cone constraints:

                   x(1)^2 >= sum(x([2 4]).^2),   x(1) >= 0,

                   x(9)^2 >= sum(x([3 5 7]).^2), x(9) >= 0,

                  are specified as

                  model.cones(1).index = [1 2 4];

                  model.cones(2).index = [9 3 5 7];

 

    model.quadcon: The quadratic constraints. A struct array. When present,

                   each element in the array defines a single quadratic

                   constraint:

                         x'*Qc*x + q'*x <= beta.

                   The Qc matrix must be a square matrix whose row and

                   column counts are equal to the number of columns of A.

                   Qc must be a sparse matrix. It is stored in

                   model.quadcon.Qc. The q vector defines the linear terms

                   in the constraint. You must specify a value for q for

                   each column of A. This must be a dense vector. It is

                   stored in model.quadcon.q.  The scalar beta defines the

                   right-hand side of the constraint. It is stored in

                   model.quadcon.rhs. For example, the ith quadratic

                   constraint is specified as

                   model.quadcon(i).Qc = Qci;

                   model.quadcon(i).q = qi;

                   model.quadcon(i).rhs = betai;

 

    model.sos: The Special Ordered Set (SOS) constraints. A struct

               array.  When present, each element in the array defines

               a single SOS constraint.  A SOS constraint can be of

               type 1 or 2. This is specified via model.sos.type. A

               type 1 SOS constraint is a set of variables for which

               at most one variable in the set may take a value other

               than zero. A type 2 SOS constraint is an ordered set of

               variables where at most two variables in the set may

               take non-zero values. If two take non-zeros values, they

               must be contiguous in the ordered set.  The members of

               an SOS constraint are specified by placing their

               indices in model.sos.index. Optional weights associated

               with SOS members may be defined in model.sos.weight.

               For example, the two SOS constraints:

               x(1) ~= 0 or x(2) ~= 0, x(1) ~= 0 or x(3) ~= 0,

               are specified as

               model.sos(1).index = [1 2];

               model.sos(1).type = 1;

               model.sos(1).weight = [1 2];

               model.sos(2).index = [1 3];

               model.sos(2).type = 1;

               model.sos(2).weight = [1 2];

 

    model.pwlobj: The piecewise-linear objective functions. A struct array.

                  When present, each element in the array defines a

                  piecewise-linear objective function of a single variable.

                  The index of the variable whose objective function is being

                  defined is stored in model.pwlobj.var. The x values for the

                  points that define the piecewise-linear function are stored

                  in model.pwlobj.x. The values in the x vector must be in

                  non-decreasing order. The y values for the points that define

                  the piecewise-linear function are stored in model.pwlobj.y.

 

    model.start: The MIP start vector. The MIP solver will attempt to

                 build an initial solution from this vector. When present,

                 you must specify a start value for each variable. This

                 must be a dense vector. Note that you can leave the

                 start value for a variable undefined---the MIP solver

                 will attempt to fill in values for the undefined start

                 values. This may be done by setting the start value

                 for that variable to nan. For example, to set the

                 start values for the first two variables and leave

                 the rest undefined do:

                 model.start = nan(size(model.A,2),1);

                 model.start(1:2) = [1; 0];

 

    model.varnames: The variable names. A cell array of strings. When

                    present, each element of the array defines the

                    name of a variable. You must specify a name for

                    each column of A.

 

    model.constrnames: The constraint names. A cell array of

                       strings. When present, each element of the

                       array defines the name of a constraint. You must

                       specify a name for each row of A.

 

 

   The params struct contains Gurobi parameters. A full list may be

   found on the Parameter page of the reference manual:

      http://gurobi.com/documentation/7.5/reference-manual/

   For example:

    params.outputflag = 0;          % Silence gurobi

    params.resultfile = 'test.mps'; % Write out problem to MPS file

    params.method     = 1;          % Use dual simplex method

 

   The gurobi function returns a struct result, with various results

   stored in its named components. The specific results that are available

   depend on the type of model that was solved, and the status of

   the optimization.

 

   The result struct will always contain the following field:

 

    result.status:  The status of the optimization, returned as a string.

                    The desired result is "OPTIMAL", which indicates an

                    optimal solution to the model was found. Other status

                    codes are possible (for example, if the model has no

                    feasible solution). See the Optimization Status

                    Codes page at

                    http://gurobi.com/documentation/7.5/reference-manual/

                    for a complete list.

 

   The result struct may contain the following fields:

 

    result.objval: The objective value of the computed solution.

 

    result.objbound: The best available bound on the solution. This

                     is a lower bound for minimization, and an

                     upper bound for maximization.

 

    result.runtime: The wall-clock time (in seconds) for the optimization.

 

    result.itercount: The number of simplex iterations performed.

 

    result.baritercount: The number of barrier iterations

                         performed.

 

    result.nodecount: The number of branch-and-cut nodes explored.

 

    result.x: The computed solution. This array contains one entry for each

              column of A.

 

    result.slack: The constraint slack in the current solution. This array

                  contains one entry for each row of A.

 

    result.qcslack: The quadratic constraint slack in the current solution.

                    This array contains one entry for each second-order

                    cone constraint and one entry for each quadratic

                    constraint. The slacks for the second-order cone

                    constraints appear before the slacks for the

                    quadratic constraints.

 

    result.rc: The reduced cost in the current solution. This is only

               available for continuous models. This array contains one

               entry for each column of A.

 

    result.pi: The dual values for the computed solution (also known as

               shadow prices). This array contains one entry for each row

               of A.

 

    result.qcpi: The dual values associated with the quadratic contraints.

                 Only available when the parameter qcpdual is set to 1.

                 This array contains one entry for each second-order

                 cone constraint and one entry for each quadratic

                 constraint. The dual values for the second-order cone

                 constraints appear before the dual values for the

                 quadratic constraints.

 

    result.vbasis: The variable basis status values for the computed

                   optimal basis. You generally should not concern yourself

                   with the contents of this array. If you wish to use an

                   advanced start later, copy this array to model.vbasis.

                   This array contains one entry for each column of A.

 

    result.cbasis: The constraint basis status values for the computed

                   optimal basis. If you wish to use an advanced start

                   later, copy this array to model.cbasis. This array

                   contains one entry for each column of A.

 

    result.unbdray: An unbounded ray (for unbounded linear models only).

                    This is a vector that, when added to any feasible

                    solution, yields a new solution that is also feasible

                    but improves the objective. Only available when

                    parameter infunbdinfo is set to 1.

 

    result.farkasdual: A Farkas infeasibility proof (for infeasible linear

                       models only). This is a dual unbounded vector.

                       Adding this vector to any feasible solution of

                       the dual model yields a new solution that is also

                       feasible but improves the dual objective.  Only

                       available when parameter infunbdinfo is set to 1.

 

    result.farkasproof: The magnitude of infeasibility violation in Farkas

                        infeasibility proof (for infeasible linear models

                        only). A Farkas infeasibility proof identifies a

                        new constraint, obtained by taking a linear

                        combination of the constraints in the model, that

                        can never be satisfied. (the linear combination is

                        available in the farkasdual attribute). This

                        attribute indicates the magnitude of the violation

                        of this aggregated constraint. Only available when

                        parameter infunbdinfo is set to 1.

 

   result.objbound: Best available bound on solution (lower bound

                    for minimization, upper bound for maximization).

 

   result.itercount: Number of simplex iterations performed.

 

   result.baritercount: Number of barrier iterations performed.

 

   result.nodecount: Number of branch-and-cut nodes explored.

 

  The following example shows how to form and solve the problem

 

  minimize    x +   y + 2 z

  subject to  x + 2 y + 3 z <= 4

              x +   y       >= 1

              x, y, z binary,

 

  with the Gurobi MATLAB interface.

 

  clear model;

  model.obj = [1; 1; 2];

  model.A = sparse([1 2 3; 1 1 0]);

  model.sense = ['<'; '>'];

  model.rhs = [4; 1];

  model.vtype = 'B';

 

  clear params;

  params.Presolve = 2;

  params.TimeLimit = 100;

 

  result = gurobi(model, params)

 

  disp(result.objval)

  disp(result.x)

 

  Copyright 2017, Gurobi Optimization, Inc.

 

  See also gurobi_setup, sparse, struct

[Johan Löfberg]

I simply think you are not sending the same version that you are running.

There is not version of YALMIP called 47. The latest version is 20180413

You cannot use matrix data in interp1. It has to be a sorted vector in the gridded space. Use step(:) and Cost_step(:) (but then again, why interpolate a convex quadratic cost..)

[Johan Löfberg]

So? That is not the issue. The issue is that you are trying to run a model (not the one you sent me) that simply cannot be solved using Gurobi. YALMIP is not even trying to call Gurobi, as it knows that the model is not possible to solve using Gurobi.

[lily ye]

I'm awfully sorry,should I update my Yalmip?MY is 20160930

Defenitely this model is the one I sent you ..

[lily ye]

I downloaded one GUROBI ,and YAlMIP has its own GUROBI sovler ?

[Johan Löfberg]

interp1 was not MISOCP-representable 2 years ago so there is your issue

but once again, using pwa representation of a convex quadratic objective when you are using a MISOCP solver is just weird

[Johan Löfberg]

no, yalmip hooks up to your gurobi installation

[lily ye]

I delete the options.IT still doesn't work

Exiting due to infeasibility:  90 lower bounds exceed the corresponding upper bounds.

[Johan Löfberg]

You haven't updated YALMIP. 

You still use the old version where interp1 isn't misocp representable, hence when you don't force it to use any specific solver, it uses one which is applicable, which means you will get fmincon, which is a terribly bad choice, since the model is a trivial SOCP

[lily ye]

Thank you again for your patience.So i need to update my Yalmip and delete the pwa..

[Johan Löfberg]

the code looks wrong anyway.

It should probably be (you want different approximators for every row in the data)

for t=1:time

    for i=1:3

        CP_linear(i,t)=interp1(step(i,:),Cost_step(i,:),P_convent(i,t),'graph');

        EP_linear(i,t)=interp1(step(i,:),Emission_step(i,:),P_convent(i,t),'graph');

    end

end

Also note that you cannot use the 'linear' flag, as that will lead to a general nonlinear callback formulation. You have to use 'graph' (if you know it is convex, which we know since you really should just use the quadratic terms, and thus get a representation which doesn't require integers) or 'milp' (for general nonconvex case) if you want to use a LP/MILP solver

[Johan Löfberg]

https://yalmip.github.io/command/interp1/

[lily ye]

Highly appreciate for your patience and help.

[lily ye]

Professor，

Thanks for your help .I downloaded the latest and rewrited the code. It does work now,but the result is different from the paper... 

I'm confused that why It warnings:

"fmincon stopped because it exceeded the function evaluation limit,

options.MaxFunctionEvaluations = 3000 (the default value)."

What is that?can I change the default value?

And also, In the paper it does said that

"constraint Pch*Pdch== 0  is redundant when charging and discharging efficiency are considered as in this work and the MILP model can be exactly relaxed to a linear model to reduce the computational burden."  But I cannot delete this constraint or the result could be weird enough

And I didn't use GUROBI, how do I know which solver it used automatically? 

[Johan Löfberg]

you are still not coding it as intended as obviously you are not using an socp solver (yalmip either doesn't find gurobi in your path, or your model  is still not an socp, but fmincon)

to serlect solver (do that) you set 'solver' to 'gurobi' in sdpsettings

Pch*Pdch should not be written as a bilinear equality, you should model that using binary variables since it simply says that at least one of them should be 0

[lily ye]

I didn't choose any solver.I tried to use GUROBI,but it warings not applicable

And for Pch*Pdch I should not simply set Pch*Pdch==0, but should model it this way:

dis=binvar(1,24);

ch=binvar(1,24);

[0<=P_ess_dis(t)<=P_ess_dis_max*dis(t),0<=P_ess_ch(t)<=P_ess_ch_max*ch(t),dis(t)+ch(t)==1];

But it doesn't work...

If I use Pch*Pdch==0, it would warning fmincon stopped.....

[Johan Löfberg]

The code you send makes no sense. You have both the beginning of the logic model (with dis+ch == 1 missing) and the nonconvex bilinear equalities. Removing the bilinear equalities and adding ch+dis <= 1 and the problem is solved in 0.2 with gurobi

[lily ye]

I guess there is something wrong with my quadratic cost, I should write it as a SOCP ?

sdpvar X,Y

minimize  Cost=XAX+BX+C+ DY

subject to  eX<=f 

Probably as  SOCP  :

sdpvar t

 minimize t

subject to  norm[A^(1/2)+A^(-1/2)*B]<=t ,eX<=f ,

and how about C,DY ？

[Johan Löfberg]

There is nothing wrong with your quadratic cost, and you should definitely not try to rewrite it to an SOCP, that is what you have YALMIP for.

[lily ye]

I'm sorry I didn't notice your reply.

But that's definitely not the right answer...In paper Cost should be more than 8000$.

My answer is only 292$. and some variables keep ZERO all the time. 

If I use Pch*Pdch==0,it performs a little bit better,and the answer would be 3000$(but it is not SOCP?) with in the iratation limit.

[lily ye]

Is there any constraint that I should transfer into SOCP?

[Johan Löfberg]

Once  again a statement which makes absolutely no sense

When you use a good structured SOCP model to implement this you get a cost of 292

When you use a bad model (nonconvex hence no guarantee to get a reasonable solution) you get a cost of 3000

A cost of 3000 is 10 times worse than 292. It does not perform better. fmincon fails miserably in finding a good solution

The fact that they report a cost of 8000 means that you are not implementing the model they have. You have missed crucial parts (or using wrong data etc)

If you get surprisingly many zeros, that indicates that you are missing some constraints or parts in the model. You should be able to think what kind of error would all zeros lead to, what constraint would not be satisfied, and thus figure out which part you are missing

[Johan Löfberg]

why on earth are you doing this now

    G4=G4+[P_ess_dis(t)+P_ess_ch(t)<=1];

[lily ye]

I'm sorry 

It should be G4=G4+[dis(t)+ch(t)<=1];(It still not the right answer...)

Not keeps zero,but keeps the minimum power constraint.

I will try to figure it out...

[lily ye]

I think h(x)==-1,I'm not sure wether it is right or wrong...

[Johan Löfberg]

As a start, P_load is wrong (first value is 220, not 20)

[Johan Löfberg]

...and your model for wind power is just completely off

mean_mean=[-0.0040,-0.0040,-0.0042,-0.0041,-0.0044,-0.0043,-0.0045,-0.0043,-0.0047,-0.0044,-0.0049,-0.0045,-0.0051,-0.0047,-0.0053,-0.0048,-0.0040,-0.0043,-0.0046,-0.0050,-0.0052,-0.0054,-0.0057,-0.0060]

Wind power is definitely not below zero and a factor 1000 smaller than everything else. It looks more like

mean_mean=[100 100 110 110 120 120 110 110 100 100 90 90 80 80 70 70 60 60 50 50 50 50 50 50];

Furter, you have E_max=80, which is clearly not used in the example since they store up to 200 in the figure

Fixing these obvious wrong data issues leads to problems where the cost is in the same neighbourhood. Note the same as in the paper though, but based on all the other errors you have, I assume there are many more things to fix.

[lily ye]

I'm so sorry for my silly mistake.Highly appreciate for your kindness and patience.