How to solve "Conversion to logical from optim.prob​lemdef.Opt​imizationE​quality is not possible." ?

Hello!

Hi, I have meet the problem in MATLAB using problem based optimization, my optimization variable is q_int , i have used to generate a vector "qu", but an error displays

"Conversion to logical from optim.problemdef.OptimizationEquality is not possible.

Error in Generate_position (line 4)

if q_int == 0"

a = floor(L/UAV_Speed*dt);
q_int = optimvar("q_int","LowerBound",0,"UpperBound",a);

function [qu] = Generate_position(Flight_direction,q_int,UAV_Speed,dt,L)   
direction = Flight_direction;
if q_int == 0   % here the error
    direction = 1;
end
if q_int == L
    direction = -1;
end
for k=2:T
    qu(k) = qu(k-1) + direction*UAV_Speed*dt;
    if qu(k)==0 || qu(k)==L
        direction = -direction;
    end
end

Appreciate for your help!

[Jan]

Please post a copy of the complete error message, which contains the information, which line is failing. You have this information on your screen already, but the readers have no chance to guess it.

[Maria]

@Jan already done, thank you.

[Torsten]

You didn't initialize qu(1) in your function. This will be another reason it will error.

[Maria]

@Torsten i initialize qu(1) = q_int*UAV_Speed*dt, but the same error is displayed,

the problem in the condition

if q_int == 0

[Matt J]

The error is produced by the part of the code that invokes Generate_position, which you haven't shown us. So, we can only guess as to how it is used.

If Generate_position is the objective function, it doesn't look like a reasonable one, since qu is a discontinuous function of q_int.

[Maria]

@Matt J

q_int is the optimized variable should have a value between 0 and a , but there is a condition that it should be put in qu_int because it related to all the vector qu

here is the function that used the vector qu to calculate a formula, and i gave above the code to generate qu

qu = Generate_position(T,Flight_direction,q_int,UAV_Speed,dt,L);
arg1 = argument1(N,T,BP,H,qu,qc);
function [arg1] = argument1(N,T,BP,H,qu,qc)
M = size(N,1);
arg1 = optimexpr(M,T);
for i=1:size(N,1)
    for k=1:T
        if N(i,k)==1   && BP(i,k)~=0
            
            arg1(i,k) = BP(i,k).*(H^2+(norm((qu(k))-qc(i,k))^2));
        end
    end
end
end

[Torsten]

@Maria

i initialize qu(1) = q_int*UAV_Speed*dt, but the same error is displayed,

I know. That's why I wrote: This will be another reason it will error.

[Maria]

@Torsten the problem is when i change q_int with such value , it works. That's why i think that the problem in the initialization of my optimization variable q_int , but i see in documents that we wrote variables as shown, or there is something missing?

q_int = optimvar("q_int","LowerBound",0,"UpperBound",a);

Answer 1

[Matt J]

This will solve the "Conversion to logical" error, but not all of the other issues raised by Torsten and me. In particular, Generate_position will still fail if q_int has any value other than 0 or L.

q_int = optimvar("q_int","LowerBound",0,"UpperBound",a);

fun=fcn2optimexpr(@(INPUT)computeThing(INPUT, T,Flight_direction,UAV_Speed,dt,L, N,T,BP,H,qc),...
                   q_int    );

function arg1=computeThing(q_int, T,Flight_direction,UAV_Speed,dt,L, N,T,BP,H,qc)

  qu = Generate_position(T,Flight_direction,q_int,UAV_Speed,dt,L);
  
  arg1 = argument1(N,T,BP,H,qu,qc);

end

[Torsten]

@Maria

There is an error in the code I supplied. You will have to use

a = floor(L/(UAV_Speed*dt));
q_int = optimvar('q_int','Type','Integer','LowerBound',0,'UpperBound',a);
fun = fcn2optimexpr(@(INPUT)computeThing(INPUT,T,Flight_direction,UAV_Speed,dt,L,N,T,BP,H,qc),...
                     q_int);

function arg1 = computeThing(q_int,T,Flight_direction,UAV_Speed,dt,L,N,T,BP,H,qc)
  qu = Generate_position(T,Flight_direction,q_int,UAV_Speed,dt,L,T);  
  arg1 = argument1(N,T,BP,H,qu,qc);
end  

function [qu] = Generate_position(Flight_direction,q_int,UAV_Speed,dt,L,T)   
  direction = Flight_direction;
  if q_int == 0   
    direction = 1;
  end
  if q_int == floor(L/(UAV_Speed*dt));
    direction = -1;
  end
  qu(1) = q_int*UAV_Speed*dt;
  for k=2:T
    qu(k) = qu(k-1) + direction*UAV_Speed*dt;
    if qu(k)==0 || qu(k)==L
        direction = -direction;
    end
  end
end

function [arg1] = argument1(N,T,BP,H,qu,qc)
  M = size(N,1);
  arg1 = optimexpr(M,T);
  for i=1:size(N,1)
    for k=1:T
        if N(i,k)==1   && BP(i,k)~=0
            arg1(i,k) = BP(i,k).*(H^2+(norm((qu(k))-qc(i,k))^2));
        end
    end
  end
end  

See whether Matt's code suggestion works as required.

[Maria]

@Matt J @Torsten i get this error

Error in fcn2optimexpr

Error in MyoptimizationProblem1 (line 71)

fun=fcn2optimexpr(@(INPUT)computeThing(INPUT, T,Flight_direction,UAV_Speed,dt,L, N,T,BP,H,qc),...

Caused by:

Function evaluation failed while attempting to determine output size. The function might contain an error, or might not be well-defined at the automatically-chosen point. To specify output size without function evaluation, use 'OutputSize'.

[Matt J]

Yes, we told you to expect this. "In particular, Generate_position will still fail if q_int has any value other than 0 or L."

[Maria]

@Matt J is there any solution that can replace this line, because the same error is displayed

q_int = optimvar("q_int","LowerBound",0,"UpperBound",a);

[Matt J]

We can run that line right here in the forum. As you can see, no errors:

a=5;
q_int = optimvar("q_int","LowerBound",0,"UpperBound",a),

q_int = 
  OptimizationVariable with properties:

          Name: 'q_int'
          Type: 'continuous'
    IndexNames: {{}  {}}
    LowerBound: 0
    UpperBound: 5

  See variables with show.
  See bounds with showbounds.

[Maria]

@Matt J so why it doesn't accept a condition , "if q_int ==0", optimvar will return just one value between 0 and "a" right?

[Torsten]

@Maria

q_int is integer. You didn't specify this in your definition.

It seems you didn't notice the modifications I posted made to your code.

M = size(N,1);

Are you sure about this ? The size of a scalar (N) is 1x1.

[Maria]

@Torsten i found a solution for the error , but i didn't know why q_int is 0, i ran the code many times and it doesn't change

Solving problem using intlinprog.
LP:                Optimal objective value is -75.000000.                                           


Optimal solution found.

Intlinprog stopped at the root node because the objective value is within a gap tolerance of the optimal value,
options.AbsoluteGapTolerance = 0 (the default value). The intcon variables are integer within tolerance,
options.IntegerTolerance = 1e-05 (the default value).


sol = 

struct with fields:

S: [96×90 double]
q_int: 0


fval =

75


exitflag = 

OptimalSolution

[Torsten]

Why is S a double array ? I thought it should be a binary integer array.

And do you have any reason to believe that q_int should be different from 0 ?

[Matt J]

so why it doesn't accept a condition "if q_int ==0".

There's no reason to think it wont't. Below is a simple example where it does:

q_int = optimvar("q_int","LowerBound",0,"UpperBound",1);

prob=optimproblem;
prob.Objective=fcn2optimexpr(@entropyFcn, q_int);

sol0.q_int=0.4;

sol=solve(prob,sol0)

Solving problem using fmincon.

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

sol = struct with fields:
    q_int: 0.3679

function out=entropyFcn(q_int)

 if q_int==0
    out=0; 
 else
    out= q_int*log(q_int);
 end

end

[Maria]

@Torsten S a double maybe because i multiplied by N with same dimension (N is the initial matrix ), but i think that q_int should change its value for many iterations and at the end it shows the optimal one that gives the best solution (maximum of S)

[Torsten]

But if S is a solution matrix, you somewhere define it as you define q_int with the "Type", "Integer" , lower and upper bounds options, don't you ?

Concerning q_int: I don't know why it doesn't change its value. Did you try to initialize it to an integer value different from 0 as Matt does it in his example code ?

[Maria]

@Torsten yes i did

S = optimvar("S",[M,T],"Type","integer","LowerBound",0,"UpperBound",1);

this error is displayed when i initialize q_int

Error using optim.internal.problemdef.ProblemImpl/solveImpl

The value of 'x0' is invalid. Initial point must contain values for variable 'S'.

[Matt J]

@Maria we aren't going to be able to help you troubleshoot effectively because you are only reporting fragments of what you are doing and what you see as a result. Like I've been doing, you should run your code here in the forum instead of on your local computer, so that we can see exactly what is being run and its result.

[Torsten]

The value of 'x0' is invalid. Initial point must contain values for variable 'S'.

I think - also with problem-based optimization - you cannot supply an initial value for only a part of your solution variables. So supplying a value only for q_int and not for S is not possible, I guess.

So I think you will additionally have to pass a matrix of 0's and 1's for x0.S (and maybe additional solution variables I don't know of).

[Torsten]

In function "Generate_position":

qu = zeros(1,T);

direction = Flight_direction;
if qu(1) == 0
    direction = 1;
end

If you set qu = zeros(1,T), especially qu(1) is 0.

Since you define the constraint

constr1 = qu(1) == q_int*UAV_Speed*dt ;

and qu(1) = 0, q_int must be 0.

Why did you change the code to something wrong ?

qu(1) must be set to q_int*UAV_Speed*dt, not to 0, with a variable integer value for q_int in the limits from 0 to floor(L/(UAV_Speed*dt)).

[Maria]

@Torsten

this constraint works only when i initialize the vector qu that's why i set qu = zeros(1,T).

constr1 = qu(1) == q_int*UAV_Speed*dt ;

also when i changed "generation position" the code works well, just for q_int = 0.

i didn't know what's the role of the solver here?

[Torsten]

i didn't know what's the role of the solver here?

What do you mean ?

[Maria]

@Torsten as i know the solver will find the best solution according to such variables, here my optimized variable is q_int and i put it between 0 and a , but it take the result for just one point "0" , there is not a sufficient number of iterations to check each value of q_int and extract the best one that gives the maximum of "S"

[Torsten]

as i know the solver will find the best solution according to such variables, here my optimized variable is q_int and i put it between 0 and a , but it take the result for just one point "0"

Yes, your problem definition is wrong.

The result from your setup will automatically be q_int = 0 because solving your constraint

qu(1) = q_int*UAV_Speed*dt

for q_int gives

q_int = qu(1)/(UAV_Speed*dt) = 0/((UAV_Speed*dt) = 0

So the solver doesn't have the possibility to vary q_int. You force the solver to set q_int = 0.

[Maria]

@Torsten but how can I define the constraint without initializing the vector qu? an error will display that qu is not recongnized.

[Torsten]

As I said: I don't know how to deal with errors arising from the problem-based approach.

Maybe it's time to switch the solution method if you have no experience with the problem-based approach.

[Maria]

@Torsten my problem can be solved with solver based approach?

[Torsten]

Did you read about the two approaches ?

If not, read this before proceeding:

https://de.mathworks.com/help/gads/global-decide-problem-solver.html

The problem-based formulation is transformed internally to a solver-based formulation.

Thus if your problem can be solved as a problem-based formulation, it can also be solved as a solver-based formulation.

[Maria]

@Torsten i have read about the two approaches, so i'm asking if there is function that transform from problem based to solver based, i see this function " prob2struct "but it seems that convert a description of problem based to matrices,really i'm wondering how can to deal with my problem, i have read about optimization but i'm failed to solve my problem.

[Maria]

@Torsten I want to eliminate this error before switching to another method.

Error in Generate_position (line 4)

if q_int == 0

can i add another variable to replace q_int for two conditions (q_int = 0 and q_int =floor(L/(UAV_Speed*dt))?

[Torsten]

The error message is due to the problem-based approach. I cannot help you here.

But wasn't the problem already solved by Matt's code ?

[Maria]

@Torsten no it wasn't solved therefore i changed q_int by qu(1) and it takes only 0 and it is not what i want

[Maria]

@Torsten @Matt J i want to use piecewise instead of if statement, but an error was displayed in this case if qu(k)==0 || qu(k)==L

function [qu] = Generate_position(Flight_direction,q_int,UAV_Speed,dt,L,T,a)
direction = Flight_direction;
syms direction(q_int)
direction(q_int) = piecewise(q_int == 0,1, q_int ==a,-1);
%   if (q_int == 0)  
%       direction = 1;
%    end
%    if q_int == floor(L/(UAV_Speed*dt))
%      direction = -1;
   %end
  qu(1) = q_int*UAV_Speed*dt;
  for k=2:T
    qu(k) = qu(k-1) + direction*UAV_Speed*dt;
    %mask = qu(k)==0 || qu(k)==L;
    syms direction(qu)
direction(qu(k)) = piecewise((qu(k) == 0)|| (qu(k)==L),-direction, direction); %here error

%     if qu(k)==0 || qu(k)==L
%         direction = -direction;
    %end
  end
end

[Matt J]

But wasn't the problem already solved by Matt's code ?...no it wasn't solved therefore i changed q_int by qu(1) and it takes only 0 and it is not what i want

I have reverted to your original definition of Generate_position, fixed some errors, and applied my suggestion to use fcn2optimexpr(). As you can see below, it produces output for qu successfully.

However, as I pointed out from the very beginning, it is also discontinuous, as you can see from when we evaluate with q_int=L-1e-12 instead of q_int=L. This small change in q_int produces a very different output. You need to make "direction" a continuous and differentiable function of q_int.

lambda = 1/10;  
T = 3600;   % Mission Time
dt = 1;     % controller time step
L = 1000;  
Vmin = 50; 
Vmax = 100; 
UAV_Speed = 20;  
%qu_start = 0;    
Flight_direction = 1;   
H = 80;      
G0 = -50;  
BW = 20e6;    
Pmax = 35;   
N0 = -130;   
C = 1e3 ;     
f = 2e9;    
alpha = 6e6;    
Beta = 10e-28;   
E_prop = 388944;   
E_max = 400e3;    
Q = 6;        

K = 10e8;

a = floor(L/UAV_Speed*dt);
q_int = optimvar("q_int","LowerBound",0,"UpperBound",a);

fcn=@(INPUT)Generate_position(Flight_direction,INPUT,UAV_Speed,dt,L,T);
generatePosition=fcn2optimexpr(fcn,q_int);

s.q_int=0;
qu=evaluate(generatePosition,s); qu(1:5)

ans = 1×5

0 20 40 60 80

s.q_int=L;
qu=evaluate(generatePosition,s); qu(1:5)

ans = 1×5

0 -20 -40 -60 -80

s.q_int=L-1e-12;
qu=evaluate(generatePosition,s); qu(1:5)

ans = 1×5

0 20 40 60 80

function [qu] = Generate_position(Flight_direction,q_int,UAV_Speed,dt,L,T)   

    direction = Flight_direction;
    if q_int == 0   % here the error
        direction = 1;
    end
    if q_int == L
        direction = -1;
    end

    qu = zeros(1,T); %<----Mattt J inserted
    for k=2:T
        qu(k) = qu(k-1) + direction*UAV_Speed*dt;
        if qu(k)==0 || qu(k)==L
            direction = -direction;
        end
    end

end

[Matt J]

i want to use piecewise instead of if statement,

Using piecewise() won't accomplish anything different from an if statement, even if you were to correct those errors. The Optimization Toolbox cannot optimize piecewise discontinuous functions.

[Torsten]

@Matt J

However, as I pointed out from the very beginning, it is also discontinuous, as you can see from when we evaluate with q_int=L-1e-12 instead of q_int=L. This small change in q_int produces a very different output. You need to make "direction" a continuous and differentiable function of q_int.

"q_int" is an integer optimization variable with value in between 0 and floor(L/(UAV_Speed*dt)) and "Flight_Direction" is another integer optimization variable with value +1 or -1, as far as I understood the optimization problem to be solved.

[Matt J]

"q_int" is an integer optimization variable with value in between 0 and

@Torsten Well, we need Maria to confirm that. So far, that condition does not show up anywhere in her remarks or in her code. Noteably, she has deliberately set S to be integer-valued, but not q_int.

[Maria]

@Matt J @Torsten until now i want to fix Flight_Direction at the begining in 1 or -1 and optimize just q_int , after this i will take the better solution from both.

for more clarification, Flight_Direction is the decision i will take at the start of q_int, it should be fixed one time, but its value will be changed through the condition of 0 and L

[Matt J]

but its value will be changed through the condition of 0 and L

This is critical to understand, and the English here is not clear enough. Should q_int be able to take any value between 0 and L? Or, should it only be able take one of two values, 0 or L?

[Torsten]

@Matt J

q_int should be able to take values 0,1,2,3,...,floor(L/(UAV_Speed*dt)) such that qu will start at qu(1) = q_int*UAV_Speed*dt which is some position in the discrete grid on [0 L] defined as

i* UAV_Speed*dt for 0<=i<=floor(L/(UAV_Speed*dt)).

[Maria]

@Matt J q_int should be able to take any value between 0 and floor(L/(UAV_Speed*dt)), for example if L = 500 and UAV_speed = 20, so q_int will choose any value between 0 and 25.

[Maria]

@Matt J exactly that @Torsten explains

[Matt J]

@Maria Ok, so Torsten is wrong? He is saying that q_int must be an integer between 0 and floor(L/(UAV_Speed*dt)), but you are saying that non-integer values like q_int=0.1 should be fine?

[Maria]

@Matt J no q_int is integer, it should be 0,1,2,3,4,5,...floor(L/(UAV_Speed*dt)),it should take any value from this range

[Maria]

@Torsten

here the code works well according to the initialization of q_int , but i'm asking if each time should i change q_int value by myself or it will be done by the solver?

fcn=@(INPUT)Generate_position(Flight_direction,INPUT,UAV_Speed,dt,L,T);
generatePosition=fcn2optimexpr(fcn,q_int);
s.q_int=5;
qu=evaluate(generatePosition,s); qu(1:T)
100	120	140	160	180	200	220	240	260	280	300	320	340	360	380	400	420	440 ....

[Torsten]

A solution variable is a variable the solver should optimize on its own. So interventions from your side are not necessary (and contraproductive).

You can see if something has changed concerning MATLAB's interpretation of your problem if "ga" is chosen as solver instead of "intlinprog". Your problem cannot be solved by "intlinprog".

[Maria]

@Torsten sorry i didn't understand you well, you mean that i change my problem to be solved by "ga"?

[Torsten]

If you use "solve" within the problem-based optimization approach, the suitable MATLAB solver is chosen automatically by the software (depending on the "difficulty" of your problem).

As you can see from the output so far:

Solving problem using intlinprog.
LP:                Optimal objective value is -75.000000.                                           
Optimal solution found.
Intlinprog stopped at the root node because the objective value is within a gap tolerance of the optimal value,
options.AbsoluteGapTolerance = 0 (the default value). The intcon variables are integer within tolerance,
options.IntegerTolerance = 1e-05 (the default value).
sol = 
struct with fields:
S: [96×90 double]
q_int: 0
fval =
75
exitflag = 
OptimalSolution

MATLAB chose "intlinprog" as solver because your problem formulation did not meet the formulation of your "real" problem.

If

Solving problem using ga.

appears here, you've made a step forward.

[Matt J]

Attached is my rewrite of your code. It runs and produces non-trivial output. Unfortunately, you have a problem. Your problem is a MINLP with over 1 million integer constrained variables. The only algorithms MATLAB has to deal with MINLP's do not handle that many unknown integers without exceeding memory limits. So, the demo below was only possible when I reduce T to something much smaller, in this case 36. You will probably need to consider reformulation of the problem.

[sol,fval,exitflag] =optimizeIt()

Solving problem using ga.
Optimization terminated: average change in the penalty fitness value less than options.FunctionTolerance
and constraint violation is less than options.ConstraintTolerance.

sol = struct with fields:
        S: [4×36 double]
    q_int: 472

fval = 3

exitflag = 
    SolverConvergedSuccessfully

[Maria]

@Torsten here it solved but q_int doesn't considered

Solving problem using ga.
Optimization terminated: average change in the penalty fitness value less than options.FunctionTolerance
and constraint violation is less than options.ConstraintTolerance.

sol = 

  struct with fields:

    S: [16×90 double]


fval =

    10

Unrecognized field name "q_int".

Error in MyoptimizationProblem1 (line 108)
sol.q_int

[Maria]

@Torsten i did show(prob), q_int not considered

OptimizationProblem : 

	Solve for:
       S
	where:
       S integer

	maximize :
       S(1, 36) + S(2, 40) + S(3, 46) + S(4, 48) + S(5, 49) + S(6, 66) + S(7, 73) + S(8, 78) + S(9, 78) + S(10, 79) + S(11,
       84)

[Torsten]

Sorry, but I don't know what you are doing and how the problem-based interface interprets what you are doing.

[Matt J]

Your problem cannot be solved by "intlinprog".

The problem is an integer linear program in S for each fixed q_int. That might be a good strategy for reformulating the problem. Loop over all possible q and solve the linear sub-problem in S.

[Torsten]

There are nonlinear constraints ...

Maybe @Maria could include the problem formulation here.

[Maria]

@Matt J " Loop over all possible q " yes this is what i need exactly , each time S will be solved for fixed q_int and extract the result, and at final it display the best q_int that gives the max of S.

i already did that without using optimization, i have tested each value that q_int could take and calculated S, but now i need to use an optimizer to solve this.

[Matt J]

i already did that without using optimization

What you describe can't be done without optimization. Even for fixed q_int, the solution for S is a non-trivial integer linear program.

In any case, I attached a solution that jointly optimizes in S and q_int, so you should check that out. However, I don't think the joint optimization can be scaled to large T.

[Maria]

@Matt J i already see it, but S is wrong, S should have just "one" by row.

[Matt J]

S should have just "one" by row.

That's not clear. Do you mean each row S(i,:) of S should contain only a single 1 and the rest should be zero? You have no such constraint in your original code. It is easy enough for you to add, however.

[Maria]

@Matt J Do you mean each row S(i,:) of S should contain only a single 1 and the rest should be zero?Yes exactly

i have N is the initial matrix that contains also a single 1 by row, so if this "1" satisfied the two constraints so S will take "1" otherwise "0"

That's why in my code i add the condition "if N(i,j)==1"

[Torsten]

so if this "1" satisfied the two constraints

??

Where do we find this in your list of constraints ?

[Matt J]

@Maria What you say implies that an optimal solution exists such that S(i,j)=0 if N(i,j)=0. However, that does not mean it is the unique optimal solution. There can be other optimal solutions which violate this, so long as (C6) and (C7) are satisfied. You may be accustomed to seeing simpler solutions for S because linear programming tends to find sparse solutions. However, now that we are jointly optimizing over both S and q_int, a linear programming algorithm can't be used.

Incidentally, your objective function implementation is unnecessarily complicated. Since N is a binary matrix, then the following is equivalent to the single expression f=S(:).'*N(:). You should probably switch to this, as it is much more efficient. In fact, the other double loops in your functions can also be easily vectorized.

    function [f] = Objective_function(S)
    f=zeros(M,T);
        for i = 1:M
            for j = 1:T
                if N(i,j) == 1
                    f1(i,j) = S(i,j).*N(i,j);
                end
            end
        end
        f = sum(f1,"all");
    end

[Torsten]

@Matt J

N is also a binary matrix - so why do you make a distinction at all ?

[Matt J]

@Torsten I don't know what "disinction" you are referencing.

[Torsten]

I meant the if-statement in the code. But now I see it's Maria's function and not your suggestion to calculate f.

[Maria]

@Matt J in my objective function i can put just sum of S

[Torsten]

in my objective function i can put just sum of S

Not according to your objective function definition.

[Matt J]

@Maria in my objective function i can put just sum of S

I think you mean f=sum(S.*N,'all'). That's true, but f=S(:).'*N(:) might be a bit faster.

Incidentally also, if you know that the only S(i,j) that matter are the ones where N(i,j)=1, then you really have only M unknown S values, whereas you have currently formulated the problem with M*T unknown values. This greatly increases the burden on the solver. You probably should remove from the problem the S(i,j) you already know can be set to 0. That would probably make it so that ga() can be scaled up to larger T. However, I still think a loop over q_int may be more reliable.

[Maria]

@Matt J i put the condition N(i,j)=1, in order to consider just ones in N and based on it S will be built.

M is unknown , it is generated randomly so its value can be changed but T is fix.

how can i use a loop over q_int

I really wonder how to solve this problem, I have passed much times to solve it but I always failed.

[Matt J]

@Maria

i put the condition N(i,j)=1, in order to consider just ones in N and based on it S will be built.

As soon as you wrote,

S = optimvar("S",[M,T],"Type","integer","LowerBound",0,"UpperBound",1);

you told the solver that there will be M*T unknown S values that it needs to solve for. The solver knows nothing about how the sparsity of N is supposed to reduce the number of unknowns. To communicate that to the solver, you could try this,

LB=zeros(M,T);
S = optimvar("S",[M,T],"Type","integer","LowerBound",LB,"UpperBound",N);

without changing anything else.

M is unknown , it is generated randomly so its value can be changed but T is fix.

No, M is fixed and known once the optimization starts, which is all that matters.

how can i use a loop over q_int

In the mathematics that you posted, constraint (C6) can be re-arranged in the form

S(i,j)<=UglyExpression(qu)

When q_int is fixed, you can precompute the right hand side. This gives a simple upper bound on S(i,j). You can then launch an optimization problem with

S = optimvar("S",[M,T],"Type","integer","LowerBound",LB,...
                                         "UpperBound",min(N,UglyExpression(qu(q_int))));

The only other constraint on S is the single linear inequaltiy (C7). The whole thing reduces to a pretty tractable looking integer linear program.

[Maria]

@Matt J sorry for my late reply! i didn't understand the use of UglyExpression ?

In the mathematics that you posted, constraint (C6) can be re-arranged in the form

S(i,j)<=UglyExpression(qu)

[Matt J]

@Maria It is what you get when you gather everything in (C6) that doesn't depend on S to one side of the inequality.

[Maria]

@Matt J you mean that i replace the constraint by this expression, sorry i didn't get it well.

[Maria]

@Torsten i'm sorry, i didn't see your comment,

so if this "1" satisfied the two constraints

??Where do we find this in your list of constraints ?

in constraint (C6) and (C7),i multiplied by S(i,j) in order to take 1 if constraints are verified otherwise "0"

[Maria]

@Matt J i already change my code to consider S just for M. Now i'm trying to find q_int because it's still fixed to a specific value given manually.

[Matt J]

@Maria

you mean that i replace the constraint by this expression, sorry i didn't get it well.

Let's take a simpler example. Suppose s is the unknown variable in a constraint that looks like,

Q - s + F(qu) >= K*s

where K is positive. Can you see that this constraint is the same as,

s <= (Q+F(qu))/(K+1)

[Maria]

@Matt J yes i get it, so i need to write an expression as a function of qu ?

[Matt J]

Yes. It should be easy. The terms in (C6) are things you are computing already, so it should be easy to re-arrange them.

[Maria]

Hello! @Matt J

i tried to apply your suggestion, but i'm asking how can i define just one input in function that contains many other functions with many inputs, for example

function [myNlCON] = Myconstraint1(qu)  % here  i wrote just (qu) as input 

qu = Generate_position(Flight_direction,q_int,UAV_Speed,dt,L,T);
        arg1 = argument1(M,N,T,H,qu,qm);
        arg2 = argument2(arg1);
        Z = sumArgument(arg2);
        
        Transmission_Delay_cons = nonlcons(Y,N,T,alpha,BP,Z);   
           Delay_exec = DelayAndProcessing(Y,N,D,T);
       myNlCON = Transmission_Delay_cons + Delay_exec ;
       
end

and i called it in my code as, here it demands arguments,should i wrote all inputs or there is a way to just write one specific

constr1 = S >= Q/(Myconstraint1(qu));

[Torsten]

Could you check

class(Myconstraint1(qu))
size(Myconstraint1(qu))
class(Q)
size(Q)

when you run your code ?

[Matt J]

@Maria i tried to apply your suggestion, but i'm asking how can i define just one input in function that contains many other functions with many inputs

Nobody said it had to have just one input. However, if you study my ealier optimizeIt.m file, you will see I avoid passing the fixed parameters as function inputs by applying the technique described here. In particular, if Myconstraint1 is nested inside some outer function, it can see (and change!) all the variables in the workspace of that function, simplifying the code to,

function outerFunction()

Flight_direction,UAV_Speed=... %The nested functions below can use these
M,N,T,H,qm,dt,L=....
Delay_exec = DelayAndProcessing(Y,N,D,T);
   

    function [myNlCON] = Myconstraint1(q_int)  %nested in outerFunction()
    
            qu = Generate_position(q_int);
            arg1 = argument1(qu);  
            arg2 = argument2(arg1);
            Z = sumArgument(arg2);
            
            Transmission_Delay_cons = nonlcons(Z);   
              
           myNlCON = Transmission_Delay_cons + Delay_exec ;
           
    end

    function qu = Generate_position(q_int) %also nested
        ...
    end

    function arg1=argument1(qu) %also nested
        ...
    end
    function arg2=argument2(arg1) %also nested
        ...
    end
    function Z=sumArgument(arg2) %also nested
        ...
    end
    function Transmission_Delay_cons = nonlcons(Z) %also nested
        ...
    end

end

[Torsten]

@Matt J

Isn't it better in general to leave the functions outside outerFunction() if they are needed somewhere else in the code ?

[Maria]

@Matt J i put all functions in outerFunction() but they still also in ohter scripts, so an error is appeared "Unrecognized function or variable " i already see your file but you are already put all functions in one script. But for me each function is in a script.

the nested functions should not exist also in other script right? only in outerFunction()?

[Matt J]

the nested functions should not exist also in other script right? only in outerFunction()?

By definition, if a function is nested inside outerFunction(), it will be in the same mfile as outerFunction.

Using nested functions was a recommendation only. You can pass all variables around as function inputs as you were doing. However, as you can see, that can get very cumbersome, and the code hard to read.

As you can see also, not all of the functions in optimizeIt.m were nested, only the ones which were specific to computing the objective function and constraints. All of the set-up functions like DelayAndProcessing were not nested, and could have been put in separate files. However, I don't know why you would want to have things like argument1() and argument2() in separate files. They seem very specific to the optimization that you are doing.

[Maria]

@Matt J i did it now, I returned to the code you provided and i have changed size of M , the problem is solved with large T, i want to verify if the q_int displayed is exactly the first point of my vector qu or not, also i want to show my constraints in order to verify the result.

i wrote show (qu) nothing was displayed.

I also want to know if there is a number of iterations was done before extracting the final result,

[Matt J]

@Maria It sounds like you have run an optimization and now have solution structure, sol with

sol.q_int=...
sol.S=...

Is that correct? If so, you would get qu by doing

qu = Generate_position(sol.q_int)

However, I expect you will not find that qu(1)=sol.q_int, but rather qu(1)=sol.q_int*UAV_Speed*dt

[Torsten]

@Matt J

You put Generate_position in OuterFunction. So it will not be found.

[Matt J]

@Torsten Well, it depends. It won't be found if you query qu from the base workspace, but it will be found if you query it inside OuterFunction. Or, OuterFunction can return a handle to Generate_position() to the base workspace, where it could be invoked freely:

function [sol,h]=outerFunction()

Flight_direction,UAV_Speed=... %The nested functions below can use these
M,N,T,H,qm,dt,L=....
Delay_exec = DelayAndProcessing(Y,N,D,T);
   
h.Generate_position =@Generate_position;

    function [myNlCON] = Myconstraint1(q_int)  %nested in outerFunction()
    
            qu = Generate_position(q_int);
            ...
           
    end

    function qu = Generate_position(q_int) %also nested
        ...
    end
   ....

end

Or, you could move Generate_position back to its own separate file and use a nested wrapper

function [sol,h]=outerFunction()

Flight_direction,UAV_Speed=... %The nested functions below can use these
M,N,T,H,qm,dt,L=....
Delay_exec = DelayAndProcessing(Y,N,D,T);

    function [myNlCON] = Myconstraint1(q_int)  %nested in outerFunction()
    
            qu = get_qu(q_int);
            ...
           
    end

    function qu = get_qu(q_int) %nested wrapper
         qu=Generate_position(Flight_direction,q_int,UAV_Speed,dt,L,T);
    end
   ....

end

[Maria]

@Torsten @Matt J qu = Generate_position(sol.q_int) works well and it displayed the vector qu where q(1) = q_int*UAV_speed*dt.

I know I asked a lot of questions, and thanks to you I managed to solve my problem, but just if you could tell me if I can add a command line to see the number of iterations performed.

[Matt J]

If you run the solver with this syntax,

[sol,fval,exitflag,report,lambda] = solve(prob) 

the variable "report" will be a structure containing various information about what the solver did, including the number of iterations.

[Maria]

@Matt J Does the number of generations refer to the number of iterations?

report = 

struct with fields:

problemtype: 'integerconstraints'
rngstate: [1×1 struct]
generations: 68
funccount: 6566
message: 'Optimization terminated: average change in the penalty fitness value less than options.FunctionTolerance↵and constraint violation is less than options.ConstraintTolerance.'
maxconstraint: 0
hybridflag: []
solver: 'ga'

[Matt J]

Yes, it does.

[Maria]

@Matt J Perfect! i really thank you for all your support ,I appreciate your taking the time.Many Thanks.

[Maria]

@Torsten I am so very thankful for your time.I am grateful for your all your support and efforts.Please accept my deepest thanks.

[Torsten]

@Maria

You are welcome.

I hope the results are as expected :-)

[Maria]

@Torsten Yes exactly! Thank you for all the help!

[Matt J]

@Maria You are quite welcome, but if you would, please Accept-click the answer to indicate that it successfully "converged".

[Maria]

hello! @Matt J

I would like to ask you if I want to change the algorithm to solve my problem to make a comparison, which algorithm could I switch on?

Thank you,

Best regards

[Matt J]

My suggestion from earlier still stands. You should try looping over q_int and solving a series of linear sub-problems.

[Maria]

@Matt J this is the method you have talking about ? the problem will be linear at each "q_int" and it can be solved with LP ?

S = optimvar("S",[M,T],"Type","integer","LowerBound",LB,...
    "UpperBound",min(N,UglyExpression(qu(q_int))));

S(i,j)<=UglyExpression(qu)

[Matt J]

Yes.

[Maria]

@Matt J I just modified this constraint, but it is still resolved by "ga".

constr1 = S >= Q./sum(fcn2optimexpr(  @Myconstraint1, q_int),2 );

[Matt J]

q_int should not be an optimvar anymore in this alternative approach. In your LP sub-problems, q_int (and therefore also qu) is known. S should be your only optimvar.

[Maria]

@Matt J q_int (and therefore also qu) is known? you mean that i put each time the value of q_int ?

if "yes", i did that but it displayed that no point satisfies the constraint

[Matt J]

Yes, that's what I mean. There's no reason to think all LP sub-problems will have non-empty constraint sets.

[Maria]

@Matt J sorry i don't understand.

[Matt J]

Your code, in outline, will look like below. Not all the LPs that you loop over will have a solution. That is why it is important to use the exitflag, as indicated, to check whether a solution was found.

for q_int=0:L
    
    i=q_int+1;
    
    S=optimvar(...)
        
    probLP=....
    
   [sol{i}, fval{i}, exitflag] = solve(probLP)
    
   if exitflag<=0 %solver failed
       fval{i}=-inf;
   end
   
   
end

[best_fval,imax]=max([fval{:}]);

solution=sol{imax};

[Torsten]

q_int goes from 0 to floor(L/(UAV_Speed*dt)) because qu(1) is defined as qu(1) = q_int*UAV_Speed*dt.

[Matt J]

The definition of qu(1) doesn't really explain why floor(L/(UAV_Speed*dt)) is the upper limit. Maybe qu is supposed to reach a maximum of L*UAV_Speed*dt.

[Torsten]

qu is supposed to reach a maximum of L and then goes back directed towards 0 by changing the variable "direction" to "-direction".

q_int as an integer variable has to satisfy 0 <= q_int <= floor(L/(UAV_Speed*dt)) for that qu(1) = q_int*UAV_Speed*dt does not exceed L.

At least that's what I took from the explanations @Maria gave so far.

[Maria]

@Matt J @Torsten it works now but i don't know how to hide this, and all statements are finished by(;)

Solving problem using intlinprog.
LP:                Optimal objective value is -53.000000.                                           


Optimal solution found.

Intlinprog stopped at the root node because the objective value is within a gap tolerance of the optimal value,
options.AbsoluteGapTolerance = 0 (the default value). The intcon variables are integer within tolerance,
options.IntegerTolerance = 1e-05 (the default value).


ans =   % how to hide this?

  1×54 cell array

  Columns 1 through 7

    {0×0 double}    {0×0 double}    {0×0 double}    {0×0 double}    {0×0 double}    {0×0 double}    {0×0 double}

  Columns 8 through 14

    {0×0 double}    {0×0 double}    {0×0 double}    {0×0 double}    {0×0 double}    {0×0 double}    {0×0 double}

  Columns 15 through 21

    {0×0 double}    {0×0 double}    {0×0 double}    {0×0 double}    {0×0 double}    {0×0 double}    {0×0 double}

  Columns 22 through 28

    {0×0 double}    {0×0 double}    {0×0 double}    {0×0 double}    {0×0 double}    {0×0 double}    {0×0 double}

  Columns 29 through 35

    {0×0 double}    {0×0 double}    {0×0 double}    {0×0 double}    {0×0 double}    {0×0 double}    {0×0 double}

  Columns 36 through 42

    {0×0 double}    {0×0 double}    {0×0 double}    {0×0 double}    {0×0 double}    {0×0 double}    {0×0 double}

  Columns 43 through 49

    {0×0 double}    {0×0 double}    {0×0 double}    {0×0 double}    {0×0 double}    {0×0 double}    {0×0 double}

  Columns 50 through 54

    {0×0 double}    {0×0 double}    {0×0 double}    {0×0 double}    {1×1 struct}

[Matt J]

Does the result agree with (or improve upon) the result that ga was giving you?

[Maria]

@Matt J yes.

[Torsten]

@Maria

And is it reasonable that for only one value of q_int, a solution exists ? At least that's how I interprete the {0x0 double}.

I would have thought there should always be a solution, but that you have to find the best one after you tried each possible value for q_int: [best_fval,imax]=max([fval{:}]);

[Matt J]

qu is supposed to reach a maximum of L and then goes back directed towards 0 by changing the variable "direction" to "-direction".

Here is the code for qu that I submitted as my answer, and which Maria accepted. If qu is taking steps of UAV_Speed*dt, then a direction change could only occur if L were landed upon exactly. For that to occur, L would need to be an exact integer multiple of UAV_Speed*dt. If L is always an exact multiple of UAV_Speed*dt, then there would be no need for the use of the floor() command in computing the upper limit of q_int.

function [qu] = Generate_position(q_int)   
    
      qu = zeros(1,T);
      qu(1) = q_int;
    
      direction = Flight_direction;
    
      for k=2:T
        qu(k) = qu(k-1) + direction;
        if qu(k)==0 || qu(k)==L
          direction = -direction;
        end
      end
    
      qu=qu*UAV_Speed*dt;
    
    end

[Maria]

@Torsten i don't understand why the optimal soltion displayed the best_fval for many times, i think that it will display all solutions and at final the best one.

and for {0×0 double} refers to the solution or what? how can i hide that?

[Torsten]

Yes, L is an exact multiple of UAV_Speed*dt. I only took the floor function to avoid floating point errors.

I don't know what code exactly you mean that Maria accepted, but I supplied this one:

function [qu] = Generate_position(Flight_direction,q_int,UAV_Speed,dt,L,T)   
    direction = Flight_direction;
    if q_int == 0   % here the error
        direction = 1;
    end
    if q_int == floor(L/(UAV_Speed*dt))
        direction = -1;
    end
    qu = zeros(1,T); %<----Mattt J inserted
    qu(1) = q_int*UAV_Speed*dt;
    for k=2:T
        qu(k) = qu(k-1) + direction*UAV_Speed*dt;
        if qu(k)==0 || qu(k)==L
            direction = -direction;
        end
    end
end

[Maria]

@Matt J i used the code of @Torsten and it works exactly that i want.

I don't have a problem with qu or q_int now.

[Torsten]

I don't have a problem with qu or q_int now.

Well, it will make a difference if you loop q_int =0:L or q_int = 0:floor(L/(UAV_Speed*dt)) or if you use Matt's or my calculation of qu.

[Matt J]

I have no sense anymore of which of the many code versions being proferred made the final cut. But since the LP approach and the ga approach are both producing the desired result, it would appear that all questions have been resolved.

[Maria]

@Matt J i used your provided code that you attached but in the part of qu i used @Torsten code.

Yes both "ga" and "LP" are producing the desired result.

Thank you very much.

[Matt J]

You're welcome. Which method is faster?

[Maria]

@Matt J LP is faster than "ga".

[Maria]

@Matt J @Torsten Hi!

hope that you are doing well.

always in the same problem, i tried to generate "ga" for a number of iterations, i ran my program but it takes too much time without given result. I'm asking if there is a mistake in my code.

this is my code :

function [sol,fval,exitflag] = outerFunction()
 %variables

 data1 = load('all_N.mat');   %load data
 all_N = data1.all_N;
 data2 = load('all_qm.mat');
 all_qm = data2.all_qm;
num = 5; %number of iteration
for iter = 1:num  
    
     qm = all_qm{iter};
     N = all_N{iter};
     
     %rest of code
     
     [sol{iter},fval{iter},exitflag,report,lambda] = solve(prob) 
end
% nested functions
end

Thanks in advance,

Regards

[Matt J]

@Maria Since this is a different question, you should probably start a new thread for it.

[Torsten]

@Maria

And of course your input files are required.

[Maria]

@Matt J @Torsten okay, thank you!