Division by an OptimizationVariable not supported.

Hi everyone,

I met a problem that in my equation, there is a part where includes a division of two optimization varaibles. For example, if i let

y=optimvar('y',[4,1,4],'Type','integer','LowerBound',0,'UpperBound',1);

then my equation would be like: (1-yi)/sum(yi)

However, matlab says division of two optimization variable is not supported by optimization toolbox, while I don't think i can change the form of my equation. So I

'm wondering if it's possible to proceed.

Thank you!

Answer 1

[Matt J]

The problem-based solver is largely limited to linear and (in R2019) quadratic programs. If you find yourself computing non-linear/non-quadratic OptimizationExpressions, it will typically mean that you have a more general optimization problem than that and must resort to another solver.

Or, you could try to reformulate your problem as a linear problem, but to know if that is possible, we would have to see the full problem, not just a single expression that occurs in it somewhere.

[BOWEN LI]

Thank you for your answer! I can attach my code here, the structure of my problem is straightforward, just the equation part may be little complicated:

%create optimization variables
y=optimvar('y',[4,1,4],'Type','integer','LowerBound',0,'UpperBound',1);
yi=[y(1,1,1),y(2,1,1),y(3,1,1),y(4,1,1);y(2,1,1),y(2,1,1),y(3,1,1),y(4,1,1);y(3,1,1),y(3,1,1),y(3,1,1),y(4,1,1);y(4,1,1),y(4,1,1),y(4,1,1),y(4,1,1)]; %create binary decision variable yit
yi(:,:,2)=[y(1,1,2),y(2,1,2),y(3,1,2),y(4,1,2);y(2,1,2),y(2,1,2),y(3,1,2),y(4,1,2);y(3,1,2),y(3,1,2),y(3,1,2),y(4,1,2);y(4,1,2),y(4,1,2),y(4,1,2),y(4,1,2)];
yi(:,:,3)=[y(1,1,3),y(2,1,3),y(3,1,3),y(4,1,3);y(2,1,3),y(2,1,3),y(3,1,3),y(4,1,3);y(3,1,3),y(3,1,3),y(3,1,3),y(4,1,3);y(4,1,3),y(4,1,3),y(4,1,3),y(4,1,3)];
yi(:,:,4)=[y(1,1,4),y(2,1,4),y(3,1,4),y(4,1,4);y(2,1,4),y(2,1,4),y(3,1,4),y(4,1,4);y(3,1,4),y(3,1,4),y(3,1,4),y(4,1,4);y(4,1,4),y(4,1,4),y(4,1,4),y(4,1,4)];
s11=optimvar('s11',[1 1 4],'Type','integer','LowerBound',0,'UpperBound',1);
s12=optimvar('s12',[1 1 4],'Type','integer','LowerBound',0,'UpperBound',1);
s13=optimvar('s13',[1 1 4],'Type','integer','LowerBound',0,'UpperBound',1);
s14=optimvar('s14',[1 1 4],'Type','integer','LowerBound',0,'UpperBound',1);
s22=optimvar('s22',[1 1 4],'Type','integer','LowerBound',0,'UpperBound',1);
s23=optimvar('s23',[1 1 4],'Type','integer','LowerBound',0,'UpperBound',1);
s24=optimvar('s24',[1 1 4],'Type','integer','LowerBound',0,'UpperBound',1);
s33=optimvar('s33',[1 1 4],'Type','integer','LowerBound',0,'UpperBound',1);
s34=optimvar('s34',[1 1 4],'Type','integer','LowerBound',0,'UpperBound',1);
s44=optimvar('s44',[1 1 4],'Type','integer','LowerBound',0,'UpperBound',1);
sij=[s11(1,1,1),s12(1,1,1),s13(1,1,1),s14(1,1,1);s12(1,1,1),s22(1,1,1),s23(1,1,1),s24(1,1,1);s13(1,1,1),s23(1,1,1),s33(1,1,1),s34(1,1,1);s14(1,1,1),s24(1,1,1),s34(1,1,1),s44(1,1,1)];
sij(:,:,2)=[s11(1,1,2),s12(1,1,2),s13(1,1,2),s14(1,1,2);s12(1,1,2),s22(1,1,2),s23(1,1,2),s24(1,1,2);s13(1,1,2),s23(1,1,2),s33(1,1,2),s34(1,1,2);s14(1,1,2),s24(1,1,2),s34(1,1,2),s44(1,1,2)];
sij(:,:,3)=[s11(1,1,3),s12(1,1,3),s13(1,1,3),s14(1,1,3);s12(1,1,3),s22(1,1,3),s23(1,1,3),s24(1,1,3);s13(1,1,3),s23(1,1,3),s33(1,1,3),s34(1,1,3);s14(1,1,3),s24(1,1,3),s34(1,1,3),s44(1,1,3)];
sij(:,:,4)=[s11(1,1,4),s12(1,1,4),s13(1,1,4),s14(1,1,4);s12(1,1,4),s22(1,1,4),s23(1,1,4),s24(1,1,4);s13(1,1,4),s23(1,1,4),s33(1,1,4),s34(1,1,4);s14(1,1,4),s24(1,1,4),s34(1,1,4),s44(1,1,4)];%sij = 1 if OD pair ij is connected in time period t
%set constants for the equations
r1=0.05;%annual demand growth rate
r2=0.1;%annual demand growth rate caused by the completion of rail routes
d=[0 20 30 35].';%di distance beween links 1-4. d1=0 since link 0-1 does not exist.
Hb0=0.25;% Bus headway in time period 0 (starting time period) 0.25hrs/bus
Hr0=0;% Train headway in time perid 0 (starting time period) 0 hrs/train since rail starts
ub=10;%cost of buses in $/hour
nc=6;%number of cars per train
uc=50;%cost of rail cars in $/car
L=500;%rail maintenance cost in $/mile
k=1000;%rail construction cost in $/mile
td=0.08;%dwell time for bus and rail in hours
vb=40;% bus operating speed
vr=60;% train operating speed
u=10; %user value of time in dollars per hour
N=4;
qijTemp=randi(10,4); %create a 4x4 OD matrix with random volumes between 0-10
qij0=tril((qijTemp-diag(diag(qijTemp))),-1)+((tril(qijTemp-diag(diag(qijTemp)),-1))).' ;%create initial OD demand qij in time period 0
%create empty optimization expressions for the equations 
qij=zeros(4,4,N);
hb=optimexpr(4);
hr=optimexpr(4);
cu=optimexpr(4);
ci=optimexpr(4);
rb=optimexpr(4);
rr=optimexpr(4);
cv=optimexpr(4);
cm=optimexpr(4);
cc=optimexpr(4);
totalcost=optimexpr(9,4);
%set constriants
for i=1:1:N
    for j=1:1:N
        for t=2:1:N
yconstr=y(i,1,t)-y(i,1,t-1)>=0; %precedence constraint on y
sconstr=sij(:,:,t-1)-sij(:,:,t-1)<=0; %precedence constraint on s
linkconstr=sum(y(i+1:j,1,t))<=(j-i)*sij(i,j,t);
        end
    end
end
%set objective functions
for t=1:N
    if t==1
    qij(:,:,t)=qij0;
    yi(:,:,t)=zeros(4,4);
    sij(:,:,t)=zeros(4,4);
    else
    qij(:,:,t)=qij(:,:,t-1)*(1+r1).^t*((1/8)*(1+sij(:,:,t)).^t);
    end
    hb(t)=2*sqrt(2*(1-yi(1,:,t))*d/vb+sum((1-yi(1,:,t))*td)*ub./(sum((1-yi(:,:,t)).*qij(:,:,t)))); %bus headway
    hr(t)=2*sqrt((2*yi(1,:,t)*d/vr+sum(yi(1,:,t)*td))*nc*uc/sum(y(:,:,t).*qij(:,:,t))); %train headway
    cu(t)=sum(qij(:,:,t).*(1-yi(:,:,t)))*(hbt/2)*(u/4)+sum(qij(:,:,t).*yi(:,:,t))*(hrt/2)*(u/4); %user wait cost
    ci(t)=((diag((1-yi(:,:,t))*qij(:,:,t)))*(d+td))*(u/vb)+(diag(yi(:,:,t)*qij(:,:,t))*(d+td))*(u+vb);%user in-vehicle cost
    rb(t)=2*(1-yi(1,:,t))*d/vb+2*sum((1-yi(1,:,t))*td);%bus round trip time
    rr(t)=2*((yi(1,:,t)*d)/vr)+sum((yi(1,:,t))*td); %rail round trip time
    cv(t)=(rbt/hbt)*ub+(rrt/hrt)*nc*uc;%vehicle operating speed 
    cm(t)=(yi(1,:,t)*d)*L;%maintenance cost
    cc(t)=(yi(1,:,t)-yi(1,:,t-1))*d*k;%construction cost
    totalcost(t)=hb(t)+hr(t)+cu(t)+ci(t)+rb(t)+rr(t)+cv(t)+cm(t)+cc(t);
end

rrproblem = optimproblem;
rrproblem.objective=sum(totalcost*ones(4,1));
rrproblem.Constraints.yconstr=yconstr;
rrproblem.Constraints.sconstr=sconstr;
rrproblem.Constraints.linkconstr=linkconstr;
[sol,fval,exitflag,output] = linprog(rrproblem);

[BOWEN LI]

The error occurs at line 65. Or which solver do you think is more capable of solving my problem?

[Matt J]

Probably ga(). It is the only solver which can handle non-linear problems with integer constraints.

[BOWEN LI]

hi，I just got an idea that if it's possible to use fcn2optimexpr to convert my non-linear equations to be used by optimization toolbox?

for example, if i want to convert hb(t)

hb(t)=2*sqrt(2*(1-yi(1,:,t))*d/vb+sum((1-yi(1,:,t))*td)*ub./(sum((1-yi(:,:,t)).*qij(:,:,t))));

I may want to use fcn2optimexpr in this way:

 hb(t)=fcn2optimexpr(@(2*sqrt(2*(1-yi(1,:,t))*d/vb+sum((1-yi(1,:,t))*td)*ub/(sum((1-yi(:,:,t)).*qij(:,:,t))))),yi);

but there's something wrong with the syntax tho.

Thank you!

[Matt J]

Here is a list of Supported Operations on Optimization Expressions and Variables. You will have to go through your operations and see which are not on the list. It appears, for example, that sqrt is not there.