Application Center - Maplesoft

# Exact solving of nonlinear optimization problems

You can switch back to the summary page by clicking here.

nonl_opt.mws

Exact  solving of nonlinear optimization problems

by Aleksas Domarkas

Vilnius University, Faculty of Mathematics and Informatics,

Naugarduko 24, Vilnius, Lithuania

aleksas@ieva.mif.vu.lt

NOTE: Program Nopt solve exactly nonlinear optimization problems.

This is Maple 7  Worksheet.

Introduction

Program Nopt solve exactly (where  possible) nonlinear optimization problems.

This program employs the  function extre ma, which is based on method of Lagrange multipliers.

Calling Sequence:

Nopt(expr, constraints);

Nopt(expr, constraints, min);

Nopt(expr, constraints, max);

Parameters:

expr        - expression whose extrema are to be found

constraints - set of constraints

Program Nopt

 > restart;

 > Nopt := proc(f, cnstr::set) option `Copyright Aleksas Domarkas, 2001,2002`; local vars, L, S, SS, k, m, K, Fmax, Fmin, sol_min, sol_MIN, sol_max, sol_MAX;     vars := indets(f, name) union indets(cnstr, name);     L := map(convert, cnstr, equality);     K := combinat[choose](L);     S := NULL;     for k to nops(K) do         if solve(K[k]) <> NULL then             extrema(f, K[k], vars, 's || k');             S := S, allvalues(s || k)         end if     end do;     SS := {};     for k to nops([S]) do         if type(S[k], set) then SS := SS union evalc(S[k])         end if     end do;     SS := remove(has, SS, I);     S := NULL;     for k to nops(SS) do         if type(SS[k], set(equation)) and         map(evalb, evalf(simplify(subs(SS[k], cnstr)))) =         {true} then S := S, value(SS[k])         end if     end do;     SS := [S];     Fmax := -infinity;     for k to nops(SS) do         if evalf(Fmax) < evalf(simplify(subs(SS[k], f))) then             Fmax := simplify(value(subs(SS[k], f)));             sol_max := SS[k]         end if     end do;     sol_MAX := sol_max;     for k to nops(SS) do         if Fmax = simplify(value(subs(SS[k], f))) and         SS[k] minus sol_max <> {} then             sol_MAX := sol_MAX, SS[k]         end if     end do;     'Fmax' = simplify(expand(Fmax)), sol_MAX;     Fmin := infinity;     for k to nops(SS) do         if evalf(simplify(subs(SS[k], f))) < evalf(Fmin) then             Fmin := simplify(subs(SS[k], f));             sol_min := SS[k]         end if     end do;     sol_MIN := sol_min;     for k to nops(SS) do         if Fmin = simplify(value(subs(SS[k], f))) and         SS[k] minus sol_min <> {} then             sol_MIN := sol_MIN, SS[k]         end if     end do;     if nargs = 3 and args[3] = min then         RETURN('F[min]' = simplify(expand(Fmin)), sol_MIN)     elif nargs = 3 and args[3] = max then         RETURN('F[max]' = simplify(expand(Fmax)), sol_MAX)     else RETURN('F[min]' = simplify(expand(Fmin)), sol_MIN,         'F[max]' = simplify(expand(Fmax)), sol_MAX)     end if end proc;

Example 1

Problem:  Minimize function

 > F1:=.08*x[1]^2-.10*x[1]*x[2]-.10*x[1]*x[3]-.10*x[1]*x[4]+.16*x[2]^2- .04*x[2]*x[3]-.04*x[2]*x[4]+.35*x[3]^2+.12*x[3]*x[4]+.35*x[4]^2;

with constraints

 > cnstr1:={0 <= x[3], 0 <= x[1], 0 <= x[2], 0 <= x[4],  1000 <= .05*x[1]+.10*x[2]+.15*x[3]+.30*x[4],  x[1]+x[2]+x[3]+x[4]<=10000};

Solving problem:

 > cnstr1:=convert(cnstr1, rational);

 > F1:=convert(F1, rational);

 > s:=Nopt(F1, cnstr1, min);

Solution:

 > s[1],seq(x[k]=subs(s[2], x[k]), k=1..4);

 > evalf(%,20);

 >

Example 2

From vec_calc  Package -- Version 7.0,   extrema.mws   by Arthur Belmonte and Philip B. Yasskin

http://www.mapleapps.com/powertools/vectorcalculus/worksheets/vec_calc7.zip

Extremize the function:

 > f:=(x,y)->x*y*exp(-x^2/2-y^2/8);

inside or on the ellipse  where

 > g:=(x,y)->x^2/4 + y^2/16;

 > Nopt(f(x,y),{g(x,y)<=1},max); Nopt(f(x,y),{g(x,y)<=1},min);

Absolute maxima occur at the interior points (1,2) and (-1,-2), and the absolute minima occur at the interior points (1,-2) and (-1,2).

 > plots[contourplot3d](f(x,y), x=-2..2, y=-4..4, axes=boxed): plots[spacecurve]([2*cos(t),4*sin(t),0],t=0..2*Pi,color=black): plots[display3d](%,%%,orientation=[-90,0]);

We replace ellipse   by     .

 > g:=(x,y)->x^2 + y^2/4;

Then

 > Nopt(f(x,y),{g(x,y)<=1},max); Nopt(f(x,y),{g(x,y)<=1},min);

 > plots[contourplot3d](f(x,y), x=-2..2, y=-4..4, axes=boxed): plots[spacecurve]([cos(t),2*sin(t),0],t=0..2*Pi,color=black): plots[display3d](%,%%,orientation=[-90,0]);

Absolute maxima occur at the boundary points ( ) and ( ) , and the absolute minima occur at the boundary points  ( )  and  ( ) .

Example 3

Problem

From Calculus and Analytic Geometry III, Fall 1997,  M.Kawski

ftp://math.la.asu.edu/pub/kawski/MAPLE/272/f_optim/optim.mws

 > z:=x^3+3*x*y^2-3*x+2;

 > mypoints:=[[1.25, -1], [1, 1.1], [-4, 3], [-1, -3], [1.25, -1]]; plot(mypoints,thickness=3,color=magenta);

We now want to find the maximum and the minimum value of z as x and y range over this polygon
(its interior, edges, and corners).

 >

Construction set of linear constraints

 > genc:=proc(A,B) local AA,BB,l,eq; geometry[point](AA,A);geometry[point](BB,B); geometry[line](l,[AA,BB],[x,y]); eq:=geometry[Equation](l); if subs(x=0,y=0,lhs(eq))>=0 then RETURN(lhs(eq)>=0) else RETURN(lhs(eq)<=0) end if; end proc;

 > n:=nops(mypoints)-1:

 > for k to n do t||k:=genc(mypoints[k],mypoints[k+1]) end do;

 > constr:={seq(t||k, k=1..n)};

 > plots[inequal](constr,x=-4..2, y=-3..3,     optionsfeasible=(color=tan),     optionsclosed=(color=green, thickness=3),     optionsexcluded=(color=yellow));

 >

Solving problem

 > c:=convert(constr,rational);

 > Nopt(z,c,min); sol_min:=%: Nopt(z,c,max);

 > sol_max:=evalf(%);

 > mypoints3d:=[seq([op(mypoints[k]),0],k=1..n+1)];

 > pol:=plots[polygonplot3d](mypoints3d, color=tan):

 > plots[contourplot3d](z, x=-4..2, y=-3..3,contours=30, axes=boxed): plottools[point](subs(sol_min[2],[x,y,0]), color=red, symbol=diamond):

 > plottools[point](subs(sol_max[2],[x,y,0]), color=red, symbol=diamond):

 > plots[display3d](pol, %%%, %%, %, orientation = [-150,45]);

 > z:='z':

 >

Example 4

From: IMSL MATH/LIBRARY, Fortran subroutines for Mathematical Applications, vol. 1 and 2,

Chapter 8: Optimization (math.pdf  4.9 M  file)

 > F4:=100*(x[2]-x[1]^2)^2+(1-x[1])^2;

 > cnstr4:={x[1]>=-2,x[1]<=0.5,x[2]>=-1,x[2]<=2};

 > cnstr4:=convert(cnstr4,rational);

 > Nopt(F4,cnstr4);

 >

Example 5

 > F5:=x^2-2*y^2+4*x*y-6*x-1;

 > cnstr5:={x>=0,y>=0,x+y<=3};

 > Nopt(F5,cnstr5);

 >

Example 6

 > F6:=x^2+y^2-12*x+16*y;

 > cnstr6:={x^2+y^2<=144, y<=abs(x)};

 > Nopt(F6,cnstr6);

 > #minimize(F6,location);

Note: If  expression  or constraints has   abs  function   Nopt  program may give false results.

 >

Example 7

 > F7:=x+y+z;cnstr7:={x^2+y^2<=z,z<=1};

 > Nopt(F7,cnstr7);

 >

Example 8

 > F8:=x+y;

 > cnstr8:={x+2*y>=2,x-y<=2,4*x-y>=0,x<=4,y<=6.5};

 > cnstr8:=convert(cnstr8,rational);

 > Nopt(F8,cnstr8);

 >

Example 9

 > F9:=x^2+y^2;

 > cnstr9:={(x-5)^2+(y-3)^2>=9,(x-5)^2+(y-3)^2<=36,x+y>=8,x>=0,y>=0};

 > sol:=Nopt(F9,cnstr9);

 > P1:=plottools[point](subs(sol[2],[x,y]),color=black): #min

 > P2:=plottools[point](subs(sol[4],[x,y]),color=blue):  #max

 > eqs:=map(convert,cnstr9,equality):

 > P3:=plots[implicitplot](eqs,x=-1..11,y=-3..9,axes=boxed,scaling=constrained):

 > plots[display](P1,P2,P3);

 >

Example 10

 > F10:=x^6+y^6-3*x^2+6*x*y-3*y^2;

 > cnstr10:={x <= 2, y <= x, 0 <= y};

 > Nopt(F10,cnstr10);

 >

Example 11

We minimize and maximize distance  between plane  and  ellipsoid    .

 > Nopt(sqrt((x-u)^2+(y-v)^2+(z-w)^2),{x^2/96+y^2+z^2-1=0,3*u+4*v+12*w-288=0});

 >

Example 12

 > F12:=8*x^2+y^2-(x^2+y^2)^2;

 > cnstr12:={x^2+y^2<=3};

 > Nopt(F12,cnstr12);

 >

Example 13

 > F13:=5*x1-2*x2+2*x3-4*x4+x5-2*x6;

 > cnstr13:={2*x1-x2+x3-2*x4+x5+x6=1,-3*x1+x2+x4-x5+x6=2, -5*x1+x2-2*x3+x4-x6=3,seq(x||k>=0,k=1..6)};

 > Nopt(F13,cnstr13);

 >

Example 14

 > F14:=(x[1]+4*x[2]-3*x[3]+x[4])/(2*x[1]+x[2]+x[3]+3*x[4]);

 > cnstr14:={x[1]-2*x[2]+x[3]=2,2*x[1]+x[2]+x[4]=6, x[1]>=0,x[2]>=0,x[3]>=0,x[4]>=0};

 > Nopt(F14,cnstr14);

 >

Example 15

 > F15:=(3*x[1]-x[2])/(x[1]+x[2]);

 > cnstr15:={x[1]+x[2]+x[3]=5,-x[1]+3*x[2]+x[4]=7, 3*x[1]-x[2]+x[5]=11,x[1]>=0,x[2]>=0,x[3]>=0,x[4]>=0,x[5]>=0};

 > Nopt(F15,cnstr15,min);

 > Nopt(F15,cnstr15,max);

Example 16

 > F16:=(y^2-x^2)*exp(1-x^2+y^2);

 > cnstr16:={x^2+y^2<=4};

 > Nopt(F16,cnstr16);

 >

While every effort has been made to validate the solutions in this worksheet, Waterloo Maple Inc. and the contributors are not responsible for any errors contained and are not liable for any damages resulting from the use of this material.