how to build a bilinear model?

[katherine xue]

i'm a new one in learning yalmip. my question is simple ,but it really bothered me quite a few days.

i have a model like this:

sdpvar: a(1,m) b(1,m) c(1,m)

ob:  sum{ t(i)*(b(i)+c(i))

s.t:   t(i)=(d(i)+a(i))/sum(d(i)+a(i))

       others...

i've no idea about it.

i've read about the bmibnb theory part but did not understand.

could u please give me some tips? or a bilinear example?

thanks for your time

[Johan Löfberg]

You mean how to write that as bilinear constraints? I suspect your last constraint is a typo, i.e., there should be no indices in the denominators. It is not clear from your description if t is a given vector, or if it is a place-holder for the rational expression. I interpreted it as an auxiliary variable defining the ratio, otherwise just remove the definition of t.

a =sdpvar(1,m);
b =sdpvar(1,m);
c =sdpvar(1,m);
t =sdpvar(1,m);
Objective = sum(t.*(b+c));
Constraints = [t*sum(d + a) == d+a];


Note though that bmibnb is a solver for general nonconvex problems, it does not have to be bilinear. Fractional models, as you have, are however best to reformulate to models not involving divisions, which here leads to bilinear models

[katherine xue]

thanks u. sorry for my poor description. you are right with t. it's a weight variable decided by a(variable) and d(known values).

i've tried your idea with an error when running the objective setence "Attempt to reference field of non-structure array".

i checked the vectors and other parts carefully, and could not figure it out.

what kind of mistake does that mean?

[katherine xue]

whlie i changed the position of variables in objective :from t.*(b+c) to  (b+c).*t ,

now it is good!

why?

 

[Johan Löfberg]

It makes no sense that it should crash, so you must have forgotten to tell me something. Please post a minimal reproducible piece of code which leads to the crash

This returns zeros as expected

m = 3;

b =sdpvar(1,m);
c =sdpvar(1,m);
t =sdpvar(1,m);

t.*(b+c) -  (b+c).*t

[katherine xue]

Sdpvar a b c t x are all (1,m) vectors.

d, f, g, (1,m) , three known vectors.

t =(a+d)/sum(a+d) which is described in constraint : t* sum(a+d) = = (a+d),sum(t)==1

c is nolinear,

c= max( b+x.*f, g ),

objective=sum( ( b+sum(c) ).*t ) ,ok

objective=sum(t .* ( b+sum(c) ) ) , with the error:

??? Error using ==> sdpvar.times at 190

[Johan Löfberg]

Please supply an explicit piece if code that runs and reproduces rhe crash

[Johan Löfberg]

OK, managed to reproduce it. Odd bug, I will look into that.

BTW, when I say reproducible code, I mean something I can copy-paste directly without having to make any guesses. When I have to guess what you did, it can get very different (and errors that depends on possible typos etc gets lost).

This reproduces your error

m = 3;
a = sdpvar(1,m);
b = sdpvar(1,m);
t = sdpvar(1,m);
x = sdpvar(1,m);


c = max( b+x.*f,g)
% Works
objective=sum((b+sum(c)).*t)
% fails
objective=sum(t.*(b+sum(c)))

This does not (it is a very different model, although it describes the same relationships)

m = 3;
a = sdpvar(1,m);
b = sdpvar(1,m);
c = sdpvar(1,m);
t = sdpvar(1,m);
x = sdpvar(1,m);


Constraints =  [c == max( b+x.*f,g)]
% Works
objective=sum((b+sum(c)).*t)
% Works
objective=sum(t.*(b+sum(c)))

[katherine xue]

thank you!

i understand.here is my code.

m=3

a=sdpvar(1,m)

b=sdpvar(1,m)

c=sdpvar(1,m)

t=sdpvar(1,m)

x=sdpvar(1,m)

e=sdpvar(1,m)

for i=1:m

c(i)=max( b(i)+x(i)*f(1), g(1));

for j=1:n

c(i)= c(i)+max( b(i)+x(i)*f(j), g(j));

end

e(i)=b(i)+c(i)/n

end

objective=e.*t; % works

objective=t.*e; % fails

 

your given example“ it is a very different model, although it describes the same relationships”. that is the question what I have, too. I wonder if it’s ok to talk about it here or open another topic.

These are two descriptions in my program:

constraints=[t*sum(a+d) == (a+d)]; % works

t=(a+d)/sum(a+d) ; fails,

“Sigmonial scalar”

What does that means? And how the two descriptions work so different? which one should i use?

 

[Johan Löfberg]

As I said, the crash is due to a silly bug, and not related to any mistake in your code. (linear expression)*(linear expression + nonlinear operator) crashes at the moment. It has been fixed already and will be included in the next release (probably released later today). As an example, x*log(x) works, log(x)*x works, but x*(x+log(x)) fails.

A constraint c*a == b (or a*c == b, same thing) vs a *assigning* a fractional expression c = b/a can make a huge difference. With the equality constrained bilinear model, there is no risk of division by zero in the numerical solver. If you however create a model which involves a fraction such as b/a, the numerical solver might encounter a situation where the current iterate of b is zero, thus leading NAN, wnd the solver will crash most likely

Signomial term means something like x^6*y^0.2*z^-6, i.e., a generalization of a fraction b/a (i.e. b^1*a^(-1))

[Johan Löfberg]

By the way, are you sure about the logic in your loop? As it is coded now, you count add the same first element to c(i) twice. I think you either should initialize c(i) to zero, or have the inner counter start at 2..

for i=1:m
    c(i) = max( b(i)+x(i)*f(1), g(1));    
    for j=2:n        
        c(i) = c(i)+max( b(i)+x(i)*f(j), g(j));        
    end    
end

Note that you can get rid of the inner loop easily

for i=1:m
   c(i) = sum(max([b(i) + x(i)*f;g],[],1));
end

[katherine xue]

thanks Johan,you really helped me a lot!

yes, j=2:n , i realized that silly mistake some time after i posted my code. : )

 

the linear matrix variable c(i) contains so many variables ,it depends on the size of f ( i.e. n). and the cycle count m will make the finial objective more variables. which often leads to the result :

"Bilinear scalar (real, homogeneous, 684 variables)

Solver stopped prematurely.

fmincon stopped because it exceeded the iteration limit,
options.MaxIter = 1000 (the default value)."

in fact,n is for monte carlo simulation time which needs to be very very big number.

how can i handle this contradiction?

 

oh , the word i asked is "sigmonial" not "signomial" ? 

wish you a happy weeken! 

[Johan Löfberg]

sigmonial is incorrect spelling of signomial (spellt wrong in yalmip...)

I am a bit surprised you get a diagnostic code from fmincon, because YALMIP should never select fmincon to solve this problem, since the max-operators, which you use in a nonconvex fashion, will force YALMIP to model this problem as a mixed-integer nonlinear program. Hence, the only applicable solvers in YALMIP are bmibnb, bnb, and bonmin.

You will have to supply the complete reproducible code for me to check this.

BTW, the problem as you are modelling it now is inherently very (very) hard, so for large n, it will most likely not be tractable to solve it using the kind of structured global modeling approach YALMIP applies here. Perhaps you do some unnecessarily complicated modeling though, it is not uncommon that very complex models can be simplified by using a more clever model.

[katherine xue]

so you mean i should work on the mathematic problem instead using yalmip?

 

 

sdpvar :a(1,m) b(1,m) c(1,m) t(1,m) x1(1,m) x2(1,m) e(1,m)

double vector :f1(1,n) f2(1,n) g(1,n) d(1,m) h(1,m) s(1,m)

double P v1 v2 v3...v9

here's y code:

P=0;

s(1)=0;

constraints=[0<=x1(:)<=1,0<=x2(:)<=1];

for i=1:m

 c(i)=max([b(i)+x1(i)*f1+x2(i)*f2;g],[],1);

 s(i+1)=s(i)+a(i)+h(i);

 a_up=v1+h(i)*v2;

 r=x1(i)*v3+x2(i)*v4+a(i)+d(i);

 p=x1(i)*v5+x2*v6;

 p_up=(d(i)+a(i))*v7;

 r_up=(d(i)+a(i))*v8;

 P=P+p;

 constraits=[constraints,0<=a(i)<=a_up,0<=p<=p_up,0<=r<=r_up];

 e(i)=b(i)+c(i)/n;

end

constraints=[constraint,P<P_up,sum(a)==sum(h)*v9,s_lw<=s(:)<=s_up,sum(t)==1,t*sum(a+d)==(a+d)];

objective= - sum(e.*t)

 

 

[Johan Löfberg]

Once again, please send reproducible code I can run. For instance, if f1 and f2 and g are row vectors as you say, then the first expression in your loop is a vector, which doesn't make sense since you are trying to put it in a single element in c.

With reproducible code, I mean code I can copy-paste to MATLAB without a single additional piece of coding.

[Johan Löfberg]

Two other bugs

You start of with s being a double, and then start putting SDPVAR objects in this. This is not possible due to a "feature" in MATLAB object orientation framework, s will be unfortunately be NaN as you write it. A working approach is (and here using the same for e to be consistent)

s = 0;
e = [];
for i = 1:m
...
s = [s;s(i)+a(i)+h(i)];
...
e = [e;b(i)+c(i)/n];
end

(an alternative is to initialize s as an sdpvar, as you did with e, but concatenation is more logical, otherwise one might think s and e are actual decision variables)

Secondly, you have a typo, constraits should be constraints

[katherine xue]

here is the code. reproduced. m and n are decreased a lot.

maybe because the reproduce ,the running results are different from my code.but it can show that the sovler is fmincon which made you a bit surprised(see above on day 2.3)

 

clear
clc;
a=sdpvar(1,3);
b=sdpvar(1,3);
c=sdpvar(1,3);
t=sdpvar(1,3);
x1=sdpvar(1,3);
x2=sdpvar(1,3);
e=sdpvar(1,3);

f1=[1119 1274 1198 924 1274];
f2=[800 400 400 800 800];
g=[0 1 2 3 4];
d=[1477 1389 1323];
h=[0.51 1.02 0.51];

P=0;
s=640;

constraints=[0<=x1(:)<=1,0<=x2(:)<=1];

for i=1:3
c(i)=sum(max([b(i)+x1(i)*f1+x2(i)*f2;g],[],1));
s=[s;s(i)+a(i)*0.8+h(i)];
 a_up=160-h(i)*0.5;
 r=x1(i)*1119+x2(i)*800-a(i)-d(i);
 p=x1(i)*0.8+x2;
 p_up=(d(i)+a(i))*0.9;
 r_up=(d(i)+a(i))*0.3;

 P=P+p;
 constraits=[constraints,0<=a(i)<=a_up,0<=p<=p_up,0<=r<=r_up];

 e(i)=b(i)+c(i)/5;
end
constraints=[constraints,P<300,sum(a)==sum(h)*0.8,128<=s(:)<=640,sum(t)==1,t*sum(a+d)==(a+d)];
objective= - sum(e.*t)
sol = solvesdp(constraints,objective);
if sol.problem == 0
a=double(a)
b=double(b)
x1=double(x1)
x2=double(x2)   
else
    display('Hmm, something went wrong!');
    end

[Johan Löfberg]

OK, the solver chosen is BNB, and then fmincon is used as local node solver.

The problem seems to be infeasible.

BTW, you should exploit the fact that sum(a)==sum(h)*.8. If you use this explicitly,

constraints=[constraints,P<300,sum(a)==sum(h)*0.8,128<=s(:)<=640,sum(t)==1,t*sum(h)*.8+t*sum(d)==(a+d)];

you get rid of a bilinear expression, and the constraints are linear (with binaries from the modelling of max). The "only" remaining complication is the bilinear objective. Making this replacement and solving the feasibility problem reveals, once again, that the problem is infeasible (I trust the solution regarding feasibility more from a good MILP solver than BNB)