Break this file up into two pieces.  The first part should be saved under
the file name fstar.m.  The script can be saved under any name ending
with .m...

_______________________________

function y=fstar(x)
%This function returns the quantity to be minimized
%by the function fmins in the main script.  We solve
%for the volume of the can of radius x(1) and height
%x(2).  We also introduce the slack variables x(3)
%and x(4) to convert the inequality constraints on
%the price and the height into equality constraints.
%We return the negative of the volume plus a penalty
%function dealing with the equality constraints and
%parameterized by the positive number P.  This variable
%is passed to the function via the global command for
%convenience.  Thus:
global P
f=-pi*x(1)^2*x(2);
g(1)=0.1*x(2)+4*pi*0.5*x(1)+(2*pi*x(1)^2+2*pi*x(1)*x(2))*6-.25+x(3)^2;
g(2)=x(2)-0.12+x(4)^2;
y=f+P*g*g';


_________________________________

clear
echo on
format compact
%In this problem we examine constrained optimization
%using the penalty function approach.  We are asked
%to design a can which is cylindrical in shape.  You are 
%given the following costs:

%cost of side seams = $0.10/meter
%cost of end seams = $0.50/meter
%cost of sheet metal = $6.00/sq. meter

%The fabrication machine can make a maximum height can
%of 12cm. What is the largest volume can that you can
%make for $0.25?  We can use slack variables to create
%equality constraints, and use penalty functions to make
%it an unconstrained optimization problem.
pause
%And now for the solution!
%This script maximizes the volume of a can subject to
%price constraints, and the requirement that the height
%of the can must be less than 0.12 meters.  The program
%uses the penalty function approach, in which both of
%the inequality constraints are turned into equality
%constraints via slack variables, and then the equalities
%are weighted by the penalty variable P.  This algorithm
%starts with some reasonable value of P, and then
%increases it until the solution no longer changes.
%The objective function is provided in the attached
%program fstar.m.
pause
%OK, we now begin:
global P
%This is our initial guess:
x=ones(1,4)'/10;
%and our initial value of P:
P=.01;
plot(x(1),x(2),'o')
axis([0,.2,0,.2])
xlabel('radius')
ylabel('height')
set(gca,'Fontsize',18)
hold on
xt=x+1;
OPTIONS(1)=0;
OPTIONS(2)=1e-5;
while norm(x-xt)>.00001,
 xt=x;
 x=fminsearch('fstar',x);
 P=P*2;
 plot(x(1),x(2),'o')
 drawnow
 echo off
 pause
end
echo on
pause
hold off
%The optimum values are thus:
radius=x(1)
height=x(2)
volume=pi*x(1)^2*x(2)
%And the final value of the penalty function
%parameter was:
P
%so the optimum height is right up against the upper
%constraint.
pause
%As you can see from the plot, the first value of
%P causes the solution to move from the initial guess
%to a new point close to the solution.  Increasing P
%causes it to move closer and closer to the final
%value.  Note that I had to reset the tolerance in the
%fminsearch routine in order to get a sufficiently accurate
%value for the intermediate solutions.  I have to be
%careful not to set the termination of the while loop
%to a smaller number than the tolerance in fminsearch!
pause
%It is interesting to note that you could have solved
%this problem as a one parameter optimization problem!
%Because the cost is linear in the height, we can 1)
%recognize that the optimum can will cost $0.25 (rather
%than keep this as an inequality constraint) and 2) use
%this information to solve for h in terms of r.  This
%is then inserted into the volume formula, and it then
%becomes a one parameter (the radius) optimization
%problem.  Any time you can reduce a problem from a
%multidimensional problem to a one-dimensional problem,
%you have made a big gain.  This is particularly true
%since 1D problems can (at least for unimodal functions)
%be solved in a finite number of steps - guaranteed.

echo off