How to solve NLP Optimization problem ?

Hello everyone! I'm new in NLP solver, i'm wondering how to write appropriately my objective function also my variables and constraints, i have seen many videos in order to understand but i didn', so i want to share with you my optimization problem may some of volunteers could help me, and any information or suggestion will be appreciated. Thank you in advance,

[Walter Roberson]

Is this a different project than your earlier one? You have been working on a vehicle optimization problem of some kind in https://www.mathworks.com/matlabcentral/answers/1891590-how-can-i-transform-condition-into-constraint-equation?s_tid=prof_contriblnk but now you want to work on a Natural Language Processing optimization problem of some kind ???

[Maria]

@Walter Roberson hi! this is the problem, and the question i posted there is a part of this , because i didn't find solution so i posted new question sharing my problem. NLP is non linear programming solver

[Walter Roberson]

I cannot tell from the problem statement which variables are being optimized over. The X with the qu(t=1), theta, n_m, s_m notation is not clear to me.

[Walter Roberson]

One thing I noticed is (C5) talks about having to be exactly .

You spent a lot of time before talking about how qu had to bounce off of L or 0 -- so for example if it was at 7 and v is 5 and L is 10, then you would want the total distance traveled to be 5, so you would want to pass through the values 8, 9, 10, bounce to reverse direction, 9, 8 .

But the magnitude calculation given in the equation is 1 dimensional, so if you are at 7 and velocity is 5 then the only valid outcomes to consider would be 7+5 = 12 or 7-5=2 and with being limited to L = 10 then the 7+5=12 would be prohibitted -- implying that if you are a location near 0 or L and you would bounce if you traveled distance v in your current direction, then you have to reverse direction immediately, rather than traveling through to the boundary and reflecting off of it. That would be an important change to what you were discussing before.

[Maria]

The variables to be optimized is qu(t=1) and Teta.

For n_m is binary matrix already known but it is generated randomly,it will be changed each iteration ( i have my program function that generated n_m).

s_m is also a binary matrix but unknown, it is the same of n_m but depending on constraints, it means s_m is n_m satisfying constraints.

i put s_m and n_m in variables because its values will be changed not for optimization.

[Maria]

@Walter Roberson so if i will change my constraint C5 in order to be ||qu(t+1)-qu(t)||<= vmax.dt? , so in this case the velocity can be changed before acrossing the limit, I'll just add the acceleration parameter.

Is there any mistake in my answer of your question " I cannot tell from the problem statement which variables are being optimized over. The X with the qu(t=1), theta, n_m, s_m notation is not clear to me"

Please don't hesitate to tell me if something is wrong and can't be solved ,because any information is important to me.

Thank you very much for your support.

[Walter Roberson]

I do not know if there is a mistaken in that response, as I am left too confused to understand it.

The variables to be optimized is qu(t=1) and Teta.

Okay, that would be two variables. And you could constrain the qu(t=1) by using a simple lower and upper bound.

s_m is also a binary matrix but unknown

Is it input or output or temporary variable ???

it means s_m is n_m satisfying constraints

Most often constraints under-specify. For example you might say that a particular 2 x 3 matrix had to have certain row sums, and certain column sums, but that would give you (number of rows) plus (number of columns) constraints whereas you have (number of rows) times (number of columns) values in the matrix, 2+3 constraints for 2*3 entries... so the array would not be fully constrained.

[Walter Roberson]

It looks to me as if your constraints (C1) to (C7) are all given as part of whatever problem you are trying to solve.

I still do not understand what is being maximized. max X ? But I do not see any X in your equations. And how do you maximize an array such as n_m ???

I am also baffled about what all this is for . What is being computed ??

What are your inputs? What is the meaning of each of them?

[Torsten]

@Walter Roberson

X = sum_m sum_t s(m,t )* n(m,t)

and s and n are binary matrices.

But the background of the problem would of course be helpful to understand what's going on.

[Maria]

@Walter Roberson @Torsten first of all thank you for your interest,

My objective function is to maximize X which is the sum of each element of each element n_m of matrix N multplied by each element s_m of matrix S. so N and S both are binary matrix having same dimension, X is just a notation of the sum.

so N is already known, i did a function that can generate N depending on such parameters

function [N] = MyMatrix(lambda,T,L,S1,S2,S3,S4,S5) %my variables for N

end

Matrix S is unknown and should be created depending in constraints.

qu(t=1) or i can noted qu_start is the variable to be optimized,it refers to the intial position of vehicle which is extracted from vector qu[1 T], at first qu_start can be initialzed to 0. qu_start=0.

C3 is the constraint of the boundary, qu_start is limited by lb=0 and ub=L.

C4 is the direction of vehicle should be optimized just for the first departure following qu_start.

C5 the velocity constraint, it is constant.

C6 all variables here are known except "S", this constraint is nonlinear, so i should write a function

function [c ceq] = nonlinear_constraint(Q,dm_ser,Bm,N0,Pm,G0,qm,alpha,H,qu_start,T,V,dt)
%this calculation is done for each element of the matrix N (N has just one '1' per row)
end

for dt

dt = 1;     % controller time step
t_vec = 0:dt:dt*T;   % time vector

here S will get 1 if this constraintif verified , 0 otherwise.

the same for C7 ,

i have E_fly and Emax are known, E-pr is refers also to the elements of N. in this constraint , each element of S verified C6 should veriy C7 in order to keep '1' else 0.

so through a number of iteration, each time the solver will show the optimal position of qu_start that can maximize X.

[Torsten]

And what's the problem you are trying to optimize behind all this ?

[Maria]

Before doing optimization, I tried to test all positions that can qu_start start from and each time qu (all vector position) will be used to calcul the distance between vehicle and other dispositive, which is mentioned in constraint C6. The problem is in position of qu because when it will be near to dispositive, they could communicate and this is considered as "agent connected" so S =1. I'm trying to explain by written because I don't know if should I formulate this in equation or what exactly. Thank you for your patience.

[Maria]

@Torsten i think that you didn't understand the problem that's why you couldn't help me, but if you can provide to me a strcucture of problem solving using fmincon according to my problem, really i will be grateful, because i didn't know if i will put it in my main program or another file script and i will call function and variables, really i'm wondering how to do this, i have seen videos and codes here in MATLAB but i didn't understand.

any information will be appreciated.

Thank you very much.

[Torsten]

The first thing you have to realize is that you cannot use "fmincon" to solve.

You have a nonlinear optimization problem with - among other variables - an unknown binary matrix S as you claim.

So the only MATLAB code you could try to use is "ga".

Make an attempt:

• Make a list of your unknowns (with dimensions) and put them all together in one long vector
• Classify your constraints as belonging to linear equalities, linear inequalities, bound constraints, integer constraints and nonlinear constraints
• Put the constraints in the adequate fields passed to "ga" (linear equalities: Aeq, beq, linear inequalities: A, b, bound constraints: lb, ub, integer constraints:intcon and nonlinear constraints :nonlcon)

[Walter Roberson]

Please describe in words what kind of problem you are trying to solve.

You have a vehicle. The vehicle moves at a constant speed. Okay, we got that.

But we have no idea what you are talking about for "depositive" or "agent connected", or what the "controller" is in this context

We have no idea what the M or S matrices represent.

alpha means ?

Bm means ?

dm_ser means ?

dt is the time step, currently set to 1. Can this change? It is important to know if dt can ever be different than 1, because you are currently using (steps times dt) as an index for arrays

means ?

means ?

means ?

G0 means ?

H means ?

L is a length of something that the vehicle is traveling on. It seems to have to do with a line.

lambda means ?

M means ?

N0 means ?

Pm means ?

Q means ?

qm means ?

qu(t) is the current vehicle position at time t, not at index t.

qu_start is the initial vehicle position

S matrix means ?

S1, S2, S3, S4, S5 means ?

T means ? Current time maybe?

V is the constant velocity of the vehicle? Or is it the constant speed of the vehicle?

means ? Is it the same thing as V ?

You should be able to create a list of every input constant and every input variable and describe each of them

[Torsten]

Assume we set values to qu(t=1) and theta.

Doesn't it automatically follow what S is ?

So my question is: can S really be influenced independent of qu(t=1) and theta or is it determined by these two parameters ?

[Maria]

@Torsten Doesn't it automatically follow what N is ? yes qu(t=1) and theta will follow N not the inverse, it means that N will not influenced by qu(t=1) and theta.

because each iteration N will be changed but we suppose that the unmaned vehicle have an idea of what will pass during the period T, so each iteration (scenarios) will start from a pecific point.

[Walter Roberson]

lambda % is the arrival of user nodes following poisson distribution

If there are any random number generations in the code then fmincon is not suitable. fmincon is only for consistent functions.

Some of the inputs appear to be integer valued or at least discrete valued. fmincon cannot be used for discrete inputs

[Maria]

@Walter Roberson so can i change my problem to linear problem and solve it with ILP solver?

[Torsten]

According to your problem formulation, your problem is nonlinear and contains binary unknowns.

As I cannot understand your answer concerning the binary matrices N and S, I'll assume that their values cannot be directy deduced from qu(t=1) and theta.

So summarizing the answer I already gave you: you'll have to use MATLAB's "ga".

Did you in the meantime follow the steps I suggested:

• Make a list of your unknowns (with dimensions) and put them all together in one long vector
• Classify your constraints as belonging to linear equalities, linear inequalities, bound constraints, integer constraints and nonlinear constraints
• Put the constraints in the adequate fields passed to "ga" (linear equalities: Aeq, beq, linear inequalities: A, b, bound constraints: lb, ub, integer constraints:intcon and nonlinear constraints :nonlcon)

?

[Walter Roberson]

I suspect that you might want to use Problem Based Optimization

[Walter Roberson]

qu_start

I cannot tell whether qu is intended to be strictly integer? If it is then fmincon() cannot be used

Your initial flight direction, teta, appears to be +1 or -1 -- discrete values. fmincon cannot handle that. However, what you can do is optimize twice, once with an initial direction of +1 and the other time with an initial direction of -1 (that is, initial teta would not be optimized over) and take the better of the two results.

[Torsten]

@Walter Roberson

Both S and N are binary, thus discrete. I think at least one of them should be optimized.

[Walter Roberson]

@Torsten

N is randomly constructed ahead of time, not something being optimized over.

Earlier https://www.mathworks.com/matlabcentral/answers/1897285-how-to-solve-nlp-optimization-problem#comment_2572485 the OP indicated that the variables to be optimized are qu(t=1) and teta -- a list that does not include S.

That is repeated in https://www.mathworks.com/matlabcentral/answers/1897285-how-to-solve-nlp-optimization-problem#comment_2575755 -- and there it is indicated that S elements will be set if the corresponding N element was successfully serviced, so that is not an input either.

If we were to make the rule that S is only set to 1 if a node in N is successfully serviced, then I think the objective could be rewritten to sum(sum(S)) which is nnz(S) -- the total number of serviced nodes. But that would be an integer and fmincon can only optimize over continuous values.

[Walter Roberson]

WIth there being only one input value being optimized over (qu_start), assuming that we have split runs with teta constant for each run, then if two different function calls returned the same number of nodes served, then fmincon() would conclude that the function was flat and would stop optiizing. So fmincon() is simply not suitable for this situation.

[Walter Roberson]

Q %delay threshold that demanded by each user to receive their requests( it is the same for all)

Does that mean that when the drone services a request, that it takes time Q to do so? If so then each time it services the request, it can travel only (dt-Q) * v in one of the steps, rather than dt * v .

It would often make sense to say that it takes time for the drone to stop to do something; there are many fewer things that a drone can do effectively just by flying by at full speed.

It seems to me likely that the drone should not be expected to "bounce" off of the ends in a situation such as this. It does not make sense for a drone to fly in a direction that there is currently no known demand to meet -- not unless the drone cannot tell where the demand is coming from. The situation is not impossible, but is not especially likely.

One could imagine a situation in which each node sends up a randomized unique beacon when it wanted service (a different beacon each time it wanted service), and did so in a low power mode that made it impractical to triangulate the source of each beacon, so that at each step the drone has to fly over the next potential location to see if it needs service. And maybe the drone is "delivering" some kind of digital authenticator so it does not need to stop or slow down as it flies. Could be, not impossible, someone just might want to do such a thing. But in cases where the drone is able to observe the list of active requests and locations, it would make sense for the drone not to bother flying further in a direction than is needed to satisfy the furthest request in that direction.

[Maria]

@Torsten i'm trying to follow the steps you provided, but the last one when i will apply "ga" , i'm asking where exactly i put my objective function , also i read about "ga" ,it used for maximization?

Ps: i want to clarify the two last lines, i did function for the nonlinear constraint,i already have a function to calculate the d_ser(m) which is noted (DelayAndProcessing) and other function which calculate the rest of formula in constraint C6, noted DataRate, so i used them.

qu = UAV_Position(T,velocity,dir,qu1,L) ; %this is a function that calculate UAV vector position 
%---------------------
S = zeros(M,T); % this is unknown matrix , it is binary 
%========================
qu1 = optimvar("qu1","LowerBound",0,"UpperBound",L); 
%----------------------------
Aeq = zeros(1,T);
beq = zeros(1,T);
for i = 1 : T-1
    
    Aeq(i) = abs(qu(i+1)-qu(i));
    beq (i)= velocity*dt;
    
end
%--------------------------------
lb=[];
ub=[];
%----------------------------
Q = ones(M,1)*Q ;
E_max = ones(M,1)*E_max ;
K = 10e7; %very large number 

function [c1,c2] = myconstr(Q,N,Delay,S,BP,Power,G,n0,T,K)
for i = 1:M
    for j = 1:T
        
        c1 = Q(i)-(DelayAndProcessing(N,Delay,T)*S(i,j) + DataRate(V,N,BP,Power,G,n0,T)) - K(S(i,j)-1);
        c2 = sum(Energy_Pr(i,j)*S(i,j))+E_prop - E_max(i) ;
    end 
end
end
%----------------------------------------

[Walter Roberson]

The calling sequence for ga is very nearly the same as the calling sequence for fmincon()

fmincon(Function, x0, A, b, Aeq, beq, lb, ub, nonlcon, options)

compared to

ga(Function, nvars, A, b, Aeq, beq, lb, ub, nonlcon, options, intcons)

where nvars is the number of variables (not a starting point like in fmincon) and intcons is a list of indices of variables that are constrained to be integers.

ga() is a minimizer. You can use it for maximization by using

objective_to_minimize = @(x) -objective_to_maximize(x)

... just like for fmincon.

[Walter Roberson]

qu1 = optimvar("qu1","LowerBound",0,"UpperBound",L); 
%----------------------------
Aeq = zeros(1,T);

No, if you are using problem based optimization, then you do not use numeric arrays for the constraints; see https://www.mathworks.com/help/optim/ug/optim.problemdef.optimizationconstraint.html

[Torsten]

Aeq = zeros(1,T);
beq = zeros(1,T);
for i = 1 : T-1
    
    Aeq(i) = abs(qu(i+1)-qu(i));
    beq (i)= velocity*dt;
    
end

As discussed earlier, the qu(i) for i>=2 don't need to be used as optimization variables or constraints. They can simply be deduced from your values for qu(1), velocity and dt.

So these constraints are unnecessary.

Further you did not yet decide if you assume that L-qu(1) is divisible by velocity*dt or not. And how you want to handle reflection at L if it is not divisible.

[Torsten]

Yes, L-qu(1) is divisible by velocity*dt , and if it is not divisible i will define new L1, should i put this in Aeq?

What is L1 ? If it is qu(t=1), you cannot simply redefine it since it is an optimization variable. You will have to define an integer variable qint with 1<= qint <= floor(L/(velocity*dt)) and use q(t=1) = qint * velocity * dt.

[Maria]

@Torsten do you mean q(t=1) becomes 1<=qu(t=1)<=qint ?

[Maria]

@Walter Roberson I have just seen now the other comments

If we were to make the rule that S is only set to 1 if a node in N is successfully serviced, then I think the objective could be rewritten to sum(sum(S)) which is nnz(S) -- the total number of serviced nodes.yes exactly

Your initial flight direction, teta, appears to be +1 or -1 -- discrete values. fmincon cannot handle that. However, what you can do is optimize twice, once with an initial direction of +1 and the other time with an initial direction of -1 (that is, initial teta would not be optimized over) and take the better of the two results. Yes i want to fix teta each time and take the better for both scenarios.

Q is delay threshold demanded by the user, but i will consider that the drone can do the service request when it is in movement, without hovering. So i will see how many request is accepting without exceed Q.

But in cases where the drone is able to observe the list of active requests and locations, it would make sense for the drone not to bother flying further in a direction than is needed to satisfy the furthest request in that direction.

because each iteration the scenario will be changed ,so the drone will know from each point will start and each direction the optimal.

[Maria]

@Torsten i'm trying to modify my code , here i got this error in the last line

"Unable to use a value of type optim.problemdef.OptimizationExpression as an index."

Q = ones(M,1)*Q ;
E_max = ones(M,1)* E_max ;
%========================
a = floor(L/UAV_Speed*dt);
q_int = optimvar("q_int","LowerBound",0,"UpperBound",a);
S = optimvar("S",[M,T],'Type','integer',"LowerBound",0,"UpperBound",1);
%----------------------------
Aeq = qu_start;
beq = q_int*UAV_Speed*dt;
%--------------------------------
lb=[];
ub=[];
%----------------------------------------
obj = Objective_function(S);

prob = optimproblem('Objective',obj);
for i = 1:M
    for j = 1:T
const1(i,j) = alpha./(BP(i,j)*n0(i,j)*H^2+norm(qu_start+(j-1)velocity*dt*Flight_direction)-qm(i,j))<= K(S(i,j)-1);
    end 
end
% prob.Constraints.constr = const1;
% show(prob) 

[Torsten]

do you mean q(t=1) becomes 1<=qu(t=1)<=qint ?

No. Define an integer variable to be optimized named "qint" which can range between 1 and floor(L/(velocity*dt)).

This variable has to be defined as integer variable (thus intcon = 1) with lower bound 1 (lb) and upper bound floor(L/(velocity*dt)) (ub).

Then, in constraints where qu(t=1) is needed, insert qint * (velocity*dt) for this variable.

But I think the problem is too complicated for you to succeed.

Don't you have classmates that could help you ?

[Maria]

@Torsten no i don't , I really thank you for your efforts and your quick response, following your steps,i have modified the code and it seems that the problem start to be solved , just i get this error : Error using ga Fitness function must be a function handle.

intcon = 1;
a = floor(L/velocity*dt);
q_int = optimvar("q_int","LowerBound",0,"UpperBound",a);
S = optimvar("S",[M,T],'Type','integer',"LowerBound",0,"UpperBound",1);
Aeq = qu_start;
beq = q_int*velocity*dt;
lb=[];
ub=[];
%----------------------------------------
obj = Objective_function(S);

prob = optimproblem('ObjectiveSense','max');
prob.Objective = obj ;
nonlcon = Q-sum(S.*Delay_exec + Transmission_Delay_cons,2) <= K.*(sum(S,2)-1);

prob.Constraints.constr = nonlcon;
%show(prob) 
S = ga(obj,[],[],[],Aeq,beq,lb,ub,nonlcon,intcon);

[Torsten]

You mix elements of solver-based and problem-based optimization. That's not possible. Either you call "ga" directly with numerical input arguments (solver-based approach) or you formulate problem-based and call "solve".

[Walter Roberson]

... or you formulate using problem-based approach, and then call routines that converts the problem based approach into array-based struct, and then call ga() on that...

[Maria]

@Torsten @Walter Roberson after calling solve

obj = Objective_function(S);
prob = optimproblem('ObjectiveSense','max');
prob.Objective = obj ;
nonlcon = Q-sum(S.*Delay_exec + Transmission_Delay_cons,2) <= K.*(sum(S,2)-1);

prob.Constraints.constr = nonlcon;
%show(prob) 
x0.S = zeros(M,T);
[sol,fval,exitflag] = solve(prob,x0);

i get this result

In intlinprog (line 139)
In optim.problemdef.OptimizationProblem/callSolver
In optim.internal.problemdef.ProblemImpl/solveImpl
In optim.problemdef.OptimizationProblem/solve
In My_optimization_problem (line 87) 
LP:                Optimal objective value is -1.378800e+06.                                        


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).

i think that i will see the sum of "S" but it returns all elements of S in "1", also it doesn't return the varaible q_int.

[Walter Roberson]

sol.q_int

[Maria]

@Walter Roberson i already did that but it returned error, and for S returned all elements "1" , really i'm sorry for posting too many questions.

[Walter Roberson]

Please post your current version of the code

[Walter Roberson]

t_vec = 0:dt:dt*T; 

That is a vector.

for t = 1:t_vec-1

When you use a non-scalar as the bound in a : operator, the effect is the same as if you had used the first element of the array. So that code is equivalent to

for t = 1:t_vec(1)-1    

but t_vec(1) is 0, and 0-1 is -1, so the loop is equivalent to

for t = 1:-1

which will not do any iterations.

Aeq2 = qu_start;
beq2 = q_int*UAV_Speed*dt ;

qu_start = 0, so Aeq2 = 0

q_int is an optimvar(). UAV_Speed = 20, dt = 1, so beq2 = 20*q_int .

If you were to activate Aeq2, beq2 (which you do not do), then that would be the constraint that 0 = q_int*20 . Which only has a single solution, q_int = 0. Which makes the existence of q_int and the constraint pointless.

for t = 1:t_vec-1

Aeq1 = norm(qu(t+1)-qu(t));
beq1 = UAV_Speed *dt ;
end

You overwrite all of Aeq1 and beq1 in every iteration of the loop. If the loop were being run at all (which it is not because 1:-1 as described above) then you would be forming a constraint only for one step. On the other hand, you do not activate the constraint so... shrug.

Considering the qu(t+1) in the loop, I would suggest to you that the proper limit for the loop should be for t = 1 : T-1

 constraint =  E_max-(E_prop +cumsum(sum(S.*Energy_Pr,2))) <= 0;

In the run I did, TaskMatrix generated a 273 x something array, so M = 273. Energy_Pr then became a 273 x 3600 array in which each row has exactly one non-zero value. In row 1, in the test run I am doing, it is Energy_Pr(1,18) that is non-zero, and it happens to have a value of 24. So when (S.*Energy_Pr,2) is done, the output for the first row is 24*S(1,18). cumsum() of all of those similar entries gives the same thing for its leading value, so the first cumsum output is 24*S(1,18) in my test run. E_max and E_prop are constants, so for the first row, the left side of the constraint comes out as -24*S(1, 18) + 11056 . And that must be <= 0 . Subtracting the constant from both sides, -24*S(1, 18) <= -11056 . Divide both sides by -24 to get S(1,18) >= 460.6667

But...

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

So S is constraint to be either 0 or 1. And there is no value in the set {0,1} such that the value is >= 450.6667

I would suggest to you that you might have wanted to code

  
 constraint =  (E_prop +cumsum(sum(S.*Energy_Pr,2))) <= E_max;

After making those changes to loop bounds and the constraints (but not activating the Aeq* constraints), and running again, I do get through the run, and all of the S entries are 1 on output. That would suggest that you have not designed your problem correctly.

[Maria]

@Walter Roberson i think that each time nonzero elements are considered too, i tried to do nnz(S) but it doesn't work.

or i can change the non constraint

nonlcon = Q-sum(S.*Delay_exec + Transmission_Delay_cons,2) <= K.*(sum(S,2)-1);

by

3<=sum(S.*Delay_exec + Transmission_Delay_cons,2) <= Q

through this S will get 1 if the constraint verified , otherwise 0.

[Maria]

@Torsten @Walter Roberson i want to ask you, if i will linearize my constraint, can i use another solver in matlab ?

[Torsten]

As long as the solver you use does not tell you that it cannot solve your problem, the reason that the solution you get is wrong is caused by your problem formulation, not by a failure of the solver.

[Torsten]

You still mix problem-based elements

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

with solver-based elements

for t = 1:T-1
    Aeq1 = norm(qu(t+1)-qu(t));
    beq1 = UAV_Speed *dt ;
end
Aeq2 = qu_start;
beq2 = q_int*UAV_Speed*dt ;
lb = [];
ub = [];

That's not possible. Decide which optimization method you prefer.

And don't you see that you overwrite the Aeq1 and beq1 values for 1 <= t <= T-2 by those obtained for t = T-1 in the loop over t ?

Further you have to define q_int and qu_start as optimization variables. The link between them is the equality constraint

qu_start - q_int*UAV_Speed*dt == 0

[Walter Roberson]

@Torsten

qu_start = 0; according to the posted code. UAV_speed and dt are constants according to the posted code. So you have the equality constraint 0 = q_int*constant and the only solution to that is q_int = 0

If qu_start were some other constant then q_int would be constrained to be exactly equal to qu_start / (UAV_speed * dt)

[Torsten]

q_int should be an integer solution variable constrained to be between 1 and floor(L/UAV_Speed*dt)) and qu_start (=qu(t=1)) be a second solution variable defined by the equation qu_start = qu_int*UAV_Speed*dt.

[Torsten]

@Luca Ferro comment moved here:

Optimization problems have no method, property, or field named 'Equations and you tried to call it here:

These are their default properties(see docs here):

[Maria]

@Torsten yes it works, i think that i started to solve the problem but there is some issues, for example this is a function that calculate argument1,

function [arg1] = argument1(N,T,BP,H,q_int,UAV_Speed,dt,Flight_direction,qc)
arg1 = optimexpr(1,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((q_int*UAV_Speed*dt+(k-1)*UAV_Speed*dt*Flight_direction)-qc(i,k))^2));
            
        end
    end
end
end

when i did show(arg1())

it displays this result, i don't know what does the sum means?

(20000000 .* (6400 + sum(((((q_int .* 20) .* 1) + 960) - 126.6667).^2, 'all')))

[Torsten]

Why do you initialize arg1 of dimension (1,T) (arg1 = optimexpr(1,T);) and fill it as arg1(N,T) (arg1(i,k) = ...) ?

Maybe the sum comes from the call to "norm". If the sum is taken over only one element, namely ((((q_int .* 20) .* 1) + 960) - 126.6667).^2, it shouldn't matter.

[Maria]

@Torsten ir is just an error, i already modify (arg1 = optimexpr(M,T)) with M is size(N,1),

if you notice in the result shown , the "960" is given by ((k-1)*UAV_Speed*dt*Flight_direction) and my distance L =500m, so it exceed the limits.

(20000000 .* (6400 + sum(((((q_int .* 20) .* 1) + 960) - 126.6667).^2, 'all')))

but did you remember the case i have told you when the position at time slot t across the limit of distance, the direction flight should be -1 multiplied by the UAV_speed.

i'm asking how can i put a condition in order to consider it in this function(shown above (argument1)), because you told me that it doesn't need to be as constraint.

sorry for asking many questions,i appreciate all your efforts.

[Torsten]

Yes, you only need q_int and flight_direction as solution variables here.

The qu(t>1) can be computed within the function nonlcon depending on q_int, flight_direction and t.

Thus defining

qu = optimvar("qu",[1,T]);

is not necessary since you know that

qu(2) = q_int + flight_direction* dt*UAV_Speed

etc.

And if qm(t) is given externally, the norm expression can be evaluated.

[Maria]

@Torsten The qu(t>1) can be computed within the function nonlcon depending on q_int and t.

it is already computed in the function argument1,

qu(t>1) = q_int*UAV_Speed*dt+(k-1)*UAV_Speed*dt*Flight_direction.(k is the time step)

here i wrote qu(t>1) as function of qu(1) which is my optimized variable q_int.

but now i need a condition to not exceed the limits, the condition should wrote as function of t .

[Torsten]

Flight_direction is a solution variable that can change with t in your code (thus an array of length t of values +/- 1) ?

[Maria]

@Torsten so it needs a condition , because 1<q_int<a , when t > a and t<2*a ,the flight direction = -1*UAV_Speed,

also when 3*a<=t<=4*a, the flight direction = 1 until complete the total time mission T.

[Torsten]

You didn't answer my question.

Can Flightdirection change for each value of t or is it given once at the start as a fixed value in order to know if to move to the left or to the right starting from position q_int ?

[Maria]

@Torsten yes the flight direction is given once at the start as a fixed value in order to know if to move to the left or to the right starting from position q_int.

i tried to add condition in the function where q_int appears, but it doesn't work,

function [arg1] = argument1(N,T,BP,H,q_int,UAV_Speed,dt,Flight_direction,qc,L)
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
            j = 1;
            x = q_int*UAV_Speed*dt+(k-1)*UAV_Speed*dt*Flight_direction;
            while k < length(x)
                xk = x(j) ;
                arg1(i,k) = BP(i,k).*(H^2+(norm((x(j))-qc(i,k))^2));
                if xk > L
                    dx = -UAV_Speed;
                elseif xk < 0
                    dx = UAV_Speed;
                else
                    j = j+1;
                    x(j) = xk;
                    arg1(i,k) = BP(i,k).*(H^2+(norm((x(j))-qc(i,k))^2));
                end
            end
            
        end
    end
end
end

[Torsten]

I don't understand all the if and else here, but - at least according to my understanding so far - the following code should generate the qu-array depending on L, T, Flight_direction, q_int, UAV_Speed and dt:

  L = 200;
  T = 547;
  Flight_direction = -1;
  UAV_Speed = 4.0;
  dt = 1.0;

  a = floor(L/(UAV_Speed*dt));
  q_int = randi(a+1,1)-1;
  
  qu(1) = q_int*UAV_Speed*dt;
  direction = Flight_direction;
  if q_int == 0
    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

[Maria]

@Torsten thank you for the code you provided ,it works well exactly as i want, but when i put it in a function in order to call it in the formula in function "argument 1" , an error displays

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

Error in Generate_position (line 4)

if q_int == 0

knowing that i removed this line

 q_int = randi(a+1,1)-1;

because i want to use q_int given by the solver.

is there any problem?

[Torsten]

Sorry, but I always work with the solver-based optimization approach. So I cannot give you advice on how to handle this problem.

[Maria]

@Torsten many thanks for your your help and all your efforts,I'll ask another question, if anyone knows how to solve this problem.

Answer 1

[Walter Roberson]

Please leave the question open as the volunteers have put considerable effort into it.

[Maria]

@Walter Roberson okay, thank you for all your effort!