How to solve NLP Optimization problem ?

16 ビュー (過去 30 日間)

Maria 2023 年 1 月 19 日

0
リンク

この質問への直接リンク

https://jp.mathworks.com/matlabcentral/answers/1897285-how-to-solve-nlp-optimization-problem

コメント済み: Maria 2023 年 2 月 1 日

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,

63 件のコメント
61 件の古いコメントを表示61 件の古いコメントを非表示

Maria 2023 年 1 月 20 日

MATLAB Online で開く

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

Walter Roberson 2023 年 1 月 22 日

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 ?

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

Walter Roberson 2023 年 1 月 22 日

MATLAB Online で開く

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 2023 年 1 月 22 日

MATLAB Online で開く

@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
%----------------------------------------

Maria 2023 年 1 月 23 日

@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 2023 年 1 月 23 日

MATLAB Online で開く

@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 2023 年 1 月 24 日

MATLAB Online で開く

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 2023 年 1 月 27 日

編集済み: Maria 2023 年 1 月 28 日

MATLAB Online で開く

@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