Is there anyone to help me for Lagrange multipliers method?

A container with an open top is to have 10 m^3 capacity and be made of thin sheet metal. Calculate the dimensions of the box if it is to use the minimum possible amount of metal. Use Lagrange multipliers method.

f=x*y+2*x*z+2*y*z   
s.t  g=x*y*z-10=0

I wrote these codes and found this answer. Is it correct solution?

syms x y z r ll
f=x*y+2*x*z+2*y*z;
g=x*y*z-10;
L=f-r*g;
gradL=gradient(L);
[xs ys zs rs]=solve(gradL==0,[x y z r],'Real',true);
h=.01;
k=.01;
for i=1:numel(xs)
 fopt=double(subs(f,[x y z],[xs(i) ys(i) zs(i)]));
 gc=subs(g,[x y z],[xs(i)+h ys(i)+k zs(i)+ll]);
 l=double(solve(gc==0,ll));
 [a j]=min(abs(l));
 l=l(j);
 fnear=double(subs(f,[x y z],[xs(i)+h ys(i)+k zs(i)+l]));
 [xs(i) ys(i) zs(i)]
 fopt
 if fopt<fnear
 disp('min')
 elseif fopt>fnear
 disp('max')
 end
end 


ans =
 
[ (5^(1/3)*16^(2/3))/4, (5^(1/3)*16^(2/3))/4, (5^(1/3)*16^(2/3))/8]
 
 
fopt =
 
   22.1042
 
min

[Rik]

Why do you keep deleting your questions instead of responding to the comments you're getting? That is both annoying and rude.

And if it is only because you are afraid your teacher sees them, why were you posting them to a public forum in the first place?

Answer 1

[Stephan]

Hi,

your result is correct. Here is an alternative solution using lagrange multiplier:

syms l b h lambda
A = l * b + 2 * b * h + 2 * l * h; % area of needed sheet metal (thickness = 1, because we are lazy guys)
V = l * b * h - 10 == 0; % volume constraint 10 m³
L = A + lambda * lhs(V); % Langrange function
dL_dl = diff(L,l) == 0; % derivative of L with respect to l
dL_db = diff(L,b) == 0; % derivative of L with respect to b
dL_dh = diff(L,h) == 0; % derivative of L with respect to h
dL_dlambda = diff(L,lambda) == 0; % derivative of L with respect to lambda
system = [dL_dl; dL_db; dL_dh; dL_dlambda]; % build the system of equations
[l_val, b_val, h_val,lambda_val] = solve(system, [l b h lambda], 'Real', true) % solve the system of equations and display the results 
results_numeric = double([l_val, b_val, h_val,lambda_val]) % show results in a vector of data type double
V_res = l_val * b_val * h_val % calculate the resulting volume (should be 10, if correct)
A_res = l_val * b_val + 2 * b_val * h_val + 2 * l_val * h_val % calculate the needed amount (area) of sheet metal

I get the same results as you get:

l_val =
 
(5^(1/3)*16^(2/3))/4
 
 
b_val =
 
(5^(1/3)*16^(2/3))/4
 
 
h_val =
 
(5^(1/3)*16^(2/3))/8
 
 
lambda_val =
 
-(5^(2/3)*16^(1/3))/5


V_res =
 
10
 
 
A_res =

3*5^(2/3)*16^(1/3)


>> double(A_res)

ans =

   22.1042

The single results for l, b, h and lambda converted to double:

results_numeric =

    2.7144    2.7144    1.3572   -1.4736

As already shown in another question from you, you can check the results using fmincon:

% Check numeric
fun = @(x)(x(1) * x(2) + 2 * x(2) * x(3) + 2 * x(1) * x(3));
numeric_solution = fmincon(fun,[2 2.5 2],[],[],[],[],[0 0 0],[10 10 10],@nonlcon)

function [c, ceq] = nonlcon(x)
c = [];
ceq = x(1) * x(2) * x(3) - 10;
end

which gives us the same results:

Local minimum found that satisfies the constraints.

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

<stopping criteria details>


numeric_solution =

    2.7144    2.7144    1.3572

Best regards

Stephan

[madhan ravi]

+1 good use of tools!

[Stephan]

Thank you!

[Duc]

what is lhs in line 4? (L = A + lambda * lhs(V))

[Stephan]

lhs(V) means left hand site of V which is

l * b * h - 10

Answer 2

[DAI LI]

what if we have inequlity constrains, how can we deal the lagrange multiplier problem in matlab?

[Hossein Malekahmadi]

you have to turn it into equal constrains like this:

f=x*y+2*x*z+2*y*z   
s.t  g=x*y*z-10=<0

you need to add another variables to make constrains from inequal to equal like this:

f=x*y+2*x*z+2*y*z   
s.t  g=x*y*z-10+s1^2=0

then you can solve it as normal one just with extra variable that you need add the derivative of it to the code as well and the code gonna look like this:

syms l b h lambda s1
A = l * b + 2 * b * h + 2 * l * h; % area of needed sheet metal (thickness = 1, because we are lazy guys)
V = l * b * h - 10+s1^2 == 0; % volume constraint 10 m³
L = A + lambda * lhs(V); % Langrange function
dL_dl = diff(L,l) == 0; % derivative of L with respect to l
dL_db = diff(L,b) == 0; % derivative of L with respect to b
dL_dh = diff(L,h) == 0; % derivative of L with respect to h
dL_dlambda = diff(L,lambda) == 0; % derivative of L with respect to lambda
dL_ds1=diff(L,s1) == 0; % derivative for added value cause  inequality of constrains
system = [dL_dl; dL_db; dL_dh; dL_dlambda;dL_ds1]; % build the system of equations
[l_val, b_val, h_val,lambda_val,s1_val] = solve(system, [l b h lambda s1], 'Real', true); % solve the system of equations and display the results 
results_numeric = double([l_val, b_val, h_val,lambda_val,s1_val]); % show results in a vector of data type double

and the results:

results_numeric =

         0         0         0         0   -3.1623
         0         0         0         0    3.1623
    2.7144    2.7144    1.3572   -1.4736         0

[muhammad ilyas khattak khattak]

@Hossein Malekahmadi greetings hossein, I am finding comments very encouraging for myself , especially yours ... i am currently calculating lagrange baased measurement of my own problem and finding empty "results_numeric"...Though i tried to obtain much help from this link, however, still struggling to get some results.... is it possible for you to rectify and share with me possible error in my script..to work quickly, i am sharing my matlab script.... many thakns in advance once again..

[muhammad ilyas khattak khattak]

function CompResourceOptmztion26April2021()
N_Total = 40;
M_Total = 5;
thetaMultiplier = sym('thetaMultiplier',[N_Total,M_Total]);
PsiMultiplier = sym('PsiMultiplier',[N_Total,M_Total]);
OmegaMultiplier = sym('OmegaMultiplier',[N_Total,M_Total]);

F = zeros(N_Total,M_Total);
F = 20.*rand(1,M_Total) + 5;               % The RSUs VEC servers assumed computational power in GHz for each Vehicle-i 
                                           % and RSU-j interaction (total 40x5 =200)... 
for iloop = 2:N_Total
    F(iloop,:) = F(1,:);
end

solutionfsmall = zeros(N_Total,M_Total);
solutionfval = zeros(N_Total,M_Total);
Rho = 0.000013;       % 13 microsecond = 13 * 10^{-6} = 13 / 1000000 =
                      % 0.000013 seconds
                      % assign(Tau,value(2/(W+1)));       % Alternative expression.... 
                      
W = 16;
% assign(Tau,value(2/(W+1)));       % Alternative expression.... 
                     % Tau(1:5,1:40) = 2/(W+1);   % probability of transmission, with contention window = 16...
                     % Transmission Power of each vehicle is 24 dBm which is equal to 251 mW which is
                     % Tau(1:5,1:40) = 2/(W+1);   % probability of transmission, with contention window = 16...
                     % Transmission Power of each vehicle is 24 dBm which is equal to 251 mW which is
Tau = 2/(W+1);

Tc = 0.000222;  % Equation (05) alternative implementation as it is a constant value preferable not to get passed from main function...
% Tc = RTS + AIFS + Delta;  % Equation (05)      %  Tcij represent the corresponding collision period.  
% Though in further equations, it is used as Tcij, but it is Tc at every instance...not variable....
                      
P = 0.251;   % subsequenlthy equal to 0.251 W.
% The SI unit standard is Watt, convention is that dB is used for signal losses and dBm and W(watt) is used 
%for power measurement but dBM is not SI unit for power. SI unit for power
%is Watt.. folwoing are used websites for the above shared information about dBm and SI derived unit of power 
% https://en.wikipedia.org/wiki/International_System_of_Units 
% https://en.wikipedia.org/wiki/DBm
%...so conversion from 24 dBm to dB is needed.
% https://en.wikipedia.org/wiki/International_System_of_Units
solutionx = zeros(N_Total,M_Total);
x1 = zeros(N_Total,M_Total);
% x = zeros(N_Total,M_Total);

% According to
% https://www.mathworks.com/help/optim/ug/optimset.html...following is the statement...
% optimoptions is recommended instead of optimset for all solvers except fzero, fminbnd, fminsearch, and lsqnonneg.

prob = optimproblem('ObjectiveSense','maximize');
opts = optimset('Display','iter','Algorithm','sqp',...
                  'MaxFunEval',inf,'MaxIter',Inf);

% Lambda = optimvar('Lambda',40,5,'Lowerbound',0,'Upperbound',1);   % Offloading factor between 0 and 1
% fsmall = optimvar('fsmall',40,5,'Lowerbound',1,'Upperbound',5);   % % Estimated range needed computation resouce for a task by a vehicle
%                                                                   % in 1 GHz and 5 GH
% x = optimvar('x',40,5,'Lowerbound',0,'Upperbound',1);      
Lambda = sym('Lambda',[40,5]);
fsmall = sym('fsmall',[40,5]);
x = sym('x',[40,5]);
% SlctdVehiclesAftrHangarian = ObjFnc16April2021V13();   % This is obtained 'x' from function ObjFnc16April2021V13()...
ObjFncHandle = @ObjFnc17April2021V15;
solutionx = ObjFncHandle();                                                                                  
                                            
 for loopi = 1:N_Total
    for loopj = 1:M_Total
        if (solutionx(loopi,loopj) ~= 0)
            x1(loopi,loopj) = 1;
        end
    end
end
N = zeros(M_Total,1);
 for indexj = 1:M_Total
     for indexiagn = 1:N_Total
         N(indexj,1) = N(indexj,1)+x1(indexiagn,indexj);         % preparation for equation 02...
     end
 end

% THandle = @FunctionCheck04March2021V11;
% T = double(THandle());

T =[4.8086   16.3470    5.3632   29.2045   21.9415;
   22.2256   12.4061   22.5986   46.5158   20.7560;
    5.4199   26.5962   24.3252    4.0599    5.6434;
   13.1514   29.4781   11.3487   27.2695   27.6891;
   28.1107   42.6842   43.7926   38.1622   56.6938;
   46.4073   22.9845   34.7637   38.0456   24.1521;
   13.1653   24.8658   21.0495   19.7708    2.4003;
   17.0525    1.0916   13.8008   10.9176   40.1073;
    6.0684   49.4158   20.0337    8.2331   36.5436;
   13.6190   29.8625   43.2797   47.0389    5.7447;
   18.6555   33.8334   35.9882   63.7029   89.8535;
   27.0191   51.6506   18.8476   62.5034   22.0511;
   16.7527   17.1631   27.5952   30.0840   38.5862;
   26.4596   46.8389    3.3176   72.6493   40.5730;
   85.4607   26.7452   44.3335   58.2425   50.6836;
   45.1008   26.1287    5.0038   18.2506   51.3368;
   45.4442   41.1230   62.2289   46.2534   37.0860;
   72.5748   30.0466   76.4128   26.4209   16.4691;
   15.2801    5.1632   19.4780   35.1089   40.3356;
    5.8535   14.8195   73.1701   82.2782   10.6715;
   30.5916   56.2552   83.8440   18.2405   32.2989;
   35.1531   66.3858   72.5657  100.5374   46.9008;
   28.5371   26.3522    7.6884   42.6167   25.4864;
   66.0583   21.8389    1.2209    8.5753   47.5901;
   23.8421    1.0720   16.2535   21.4613   66.7389;
   37.0624   43.2530   32.1375   63.4302   72.9162;
   11.8703    1.3372   12.4381   25.4689   41.2694;
   42.2559   10.4953   31.6051   43.0764   10.1054;
   16.5871    6.8371    8.0623   28.8890    2.8208;
   10.9037   19.3624   17.3634    1.9732    6.1935;
   10.4831   39.1682   11.7113   30.8649   27.0540;
   59.9239   36.5261   17.4829   32.2847   30.6131;
   38.5933   21.4846   49.4188   46.0527   28.7820;
   78.9164   30.0680    1.6078   48.1825   55.2043;
   16.2801    2.2293    4.8801   23.7330   25.1516;
   79.5307    9.9379   24.1760    8.0106   16.8871;
   67.1426   63.4992   70.9382   17.3182   23.4616;
   35.6835    6.2871   20.8957    3.5507    9.1070;
   27.4826   37.9915   42.1662   15.6982   16.4394;
   32.9603   30.1127   79.0386   32.8151   80.9879]

Tmax = zeros(N_Total,M_Total);              % Implementation of TiMax   
                                            % Instead of (N_Total,1), TiMax is expanded to (N_Total,M_Total).....  % Tmax is now (N_Total,M_Total) matrix...
for iloop = 1:N_Total                       
   Tmax(iloop,:) = max(T(iloop,:),[],2);    % Tmax is now (N_Total,M_Total) matrix...
end

G = 2.*rand(N_Total,M_Total) + 1;   % Channel gain initially predicted to be 1 dB abnd 2 dB...
B = 10;                             % B is considered constant equal to 10 MHz  
                                    % Computation resoruce needed by each task of a vehicle... between 4 GHz and 0.5 GHz
                                      % i have to prolong the dimension of ci to get it into the equation                     
ci = zeros(N_Total,M_Total);                                      
ci(:,1) = 4.*rand(N_Total,1) + 0.5;
for jloop = 2:M_Total
    ci(:,jloop) = ci(:,1);
end

% d = 3.*rand(N_Total,1) + 1;   % 1 Mega bytes to 3 Mega bytes... 

d = zeros(N_Total,M_Total);                                      
d(:,1) = 3.*rand(N_Total,1) + 1;
for jloop = 2:M_Total
    d(:,jloop) = d(:,1);
end

func1 = FunctionCheck04March2021V11();
func1.CompTInitialHandle()
% func1.CompTsHandle;
PropTime = @ComputeTp;

Alpha1 = sym('Alpha1');
Beta2 = sym('Beta2');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [phi,Ts] = ComputeTs(indexi,indexj,Lambda,d,G,P,B)

% The wireless access technique between a vehicle and an 
% RSU is generally based on IEEE 802.11p which includes a
% contention-based medium access control (MAC) protocol,     

PHYhead = 0.000064;   % Physcial header is 64 Micro Seconds... Using SI unit of seconds for time,
                      % 64 microsecond = 64 * 10^{-6} = 64 / 1000000 =
                      % 0.000064 seconds
MAChead = 0.000043;   % 43 microsecond = 43 * 10^{-6} = 43 / 1000000 =
                      % 0.000043 seconds
ACK = 0.000103;       % 101 microsecond = 101 * 10^{-6} = 101 / 1000000 =
                      % 0.000101 seconds
Rho = 0.000013;       % 13 microsecond = 13 * 10^{-6} = 13 / 1000000 =
                      % 0.000013 seconds
                      % "Rho" is "slot" or single period time...
SIFS = 0.000032;      % 32 microsecond = 32 * 10^{-6} = 32 / 1000000 =
                      % 0.000032 seconds
AIFS = 0.000058;      % 58 microsecond = 58 * 10^{-6} = 58 / 1000000 =
                      % 0.000058 seconds
RTS = 0.000144;       % 144 microsecond = 144 * 10^{-6} = 144 / 1000000 =
                      % 0.000144 seconds
Delta = 0.000002;     % 2 microsecond = 2 * 10^{-6} = 2 / 1000000 =
                      % 0.000002 seconds
CTS = 0.000120;       % 120 microsecond = 120 * 10^{-6} = 120 / 1000000 =
                      % 0.000120 seconds
CWmin = 15*Rho;       % CWmin is 15 times slot. That is 15*Rho..
                      % That is 15*13 microseconds
                      % That is 15* 0.000013 = 0.000195 seconds
CWmax = 1023*Rho;       % CWmax is 1023 times slot. That is 1023*Rho..
                      % That is 1023*13 microseconds
                      % That is 1023* 0.000013 = 0.013299 seconds
% V = 5; % "V" is capital "V". Note that small "v" is also used in this code which is Vehicle speed.
       % "V" (capital V) in this code is used for "Maximum Backoff stage",
       % set to five (5)....
% z = 4; % "z" is "retransmission limit".
H = PHYhead + MAChead;
phi = H + SIFS + Delta + ACK+AIFS + Delta + RTS+SIFS + Delta + CTS+SIFS + Delta;
% TpHandle = @ComputeTp; 
% Ts = phi + TpHandle(Lambda,d,indexi,indexj,G,P,B);   % Equation (02)  % % Implemented inside for loops
Tc = RTS + AIFS + Delta;  % Equation (05)      %  Tcij represent the corresponding collision period.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                         % Initializing optimization parameters.....
% fsmall = 4.*rand(N_Total,M_Total) + 1;    % Estimated range needed computation resouce for a task by a vehicle
                                    % in 1 GHz and 5 GHz

x0.Alpha1 = 15;
x0.Beta2 = max(T(:,:),[],[1 2]);
x0.x(:,:) = x1(:,:);
x0.Lambda(:,:) = randi([0,1],N_Total,M_Total);
x0.fsmall(:,:) = randi([0,4],N_Total,M_Total);                      

x = x1;
% options = optimoptions('Display','iter');

for jloop = 1:M_Total
  for iloop = 1:N_Total
LambdaMatrix = Lambda(1:iloop,1:jloop);


TMovement = func1.CompMoveT(jloop);

PIdle = (1 - Tau)*N(jloop,1);  % Equation (06)

Pc = 1-(1-Tau).^N(jloop,1)-N(jloop,1).*Tau.*(1-Tau).^(N(jloop,1)-1);  % Equation (04)   % pcij denote collision probability
% pcij denote collision probability...
SuccessTrans = Tau.*N(jloop,1).*(1-Tau).^(N(jloop,1)-1);    % N(i,j) is number of ...
TMatrix = T(1:iloop,1:jloop);            % Immitating SUmmation of utility function....
xMatrix = x(1:iloop,1:jloop);
fsmallMatrix = fsmall(1:iloop,1:jloop);
LambdaMatrix = Lambda(1:iloop,1:jloop);
ciMatrix = ci(1:iloop,1);
reshape(LambdaMatrix,[iloop jloop]);     % Converting marix into column vector...
reshape(fsmallMatrix,[iloop jloop]);        % Converting marix into column vector...
reshape(xMatrix,[iloop jloop]);       % Converting marix into row vector...
reshape(TMatrix,[iloop jloop]);     % Converting marix into column vector...

constr1a = sum(x(1:N_Total,1:M_Total).*fsmall(1:N_Total,1:M_Total)) <= F(iloop,jloop);   % condition 13c RHS
constr1b = sum(x(1:N_Total,1:M_Total).*fsmall(1:N_Total,1:M_Total)) >= 0;                % condition 13c LHS
Delta = ci(iloop,1:M_Total)./fsmall(iloop,1:M_Total) + TMovement + phi + ((PIdle*Rho+Pc)*Tc)./SuccessTrans;  % Organ 2... tMovement is sum(Lk/v)....
part1 = Lambda(iloop,1:M_Total).*d(iloop,1:M_Total)/B.*log(1 + P.*G(iloop,1:M_Total));   % Psi
constr2 = sum(x(iloop,1:M_Total).*Lambda(iloop,1:M_Total).*(part1+ Delta)) <= Tmax(iloop,jloop);    % Implementation of TiMax   
                                                                                                  % Instead of (N_Total,1),TiMax is expanded to (N_Total,M_Total).....  % Tmax is now (N_Total,M_Total) matrix...
organ1 = log(1+Beta2-x(1:N_Total,1:M_Total).*Lambda(1:N_Total,1:M_Total));
ObjctFunc = sum(Alpha1.*organ1(1:N_Total,1:M_Total),'all');

LagrangianF = sum(Alpha1.*organ1(1:N_Total,1:M_Total),'all') - thetaMultiplier(iloop,jloop).*(sum(x(iloop,1:M_Total).*Lambda(iloop,1:M_Total).*(part1+ Delta)) - Tmax(iloop,jloop)) - PsiMultiplier(iloop,jloop).*(sum(x(1:N_Total,1:M_Total).*fsmall(1:N_Total,1:M_Total),'all')-F(iloop,jloop)) - OmegaMultiplier(iloop,jloop).*(sum(x(1:N_Total,1:M_Total).*fsmall(1:N_Total,1:M_Total),'all'));
dL_dthetaMultiplier = diff(LagrangianF,thetaMultiplier(iloop,jloop)) == 0;
dL_dPsiMultiplier = diff(LagrangianF,PsiMultiplier(iloop,jloop)) == 0;
dL_dOmegaMultiplier = diff(LagrangianF,OmegaMultiplier(iloop,jloop)) == 0;
dL_dLambda = diff(LagrangianF,Lambda(iloop,jloop)) == 0;
dL_dfsmall = diff(LagrangianF,fsmall(iloop,jloop)) == 0;
dL_dx = diff(LagrangianF,x(iloop,jloop)) == 0;
system1 = [dL_dLambda;dL_dfsmall;dL_dx;dL_dthetaMultiplier;dL_dPsiMultiplier;dL_dOmegaMultiplier];

[Lamba_val fsmall_val x_val thetaMultiplier_val PsiMultiplier_val OmegaMultiplier_val] = solve(system1,[Lambda(iloop,jloop) fsmall(iloop,jloop) x(iloop,jloop) thetaMultiplier(iloop,jloop) PsiMultiplier(iloop,jloop) OmegaMultiplier(iloop,jloop)],'Real',true); 
% Return only real solutions by setting 'Real' option to true.
results_numeric = double([Lamba_val,fsmall_val,x_val,thetaMultiplier_val,PsiMultiplier_val,OmegaMultiplier_val]);
    end
end
solutionfsmall
end