function ConstrainedOptimization
LW = 'linewidth'; lw = 1.6;
FS = 'fontsize'; fs = 14;
MS = 'markersize'; ms = 20;

### One-dimensional constrained optimization

By virtue of Chebfun's capabilities with global rootfinding and global optimization, some constrained optimization is possible. For example, for $x\in[0,10]$, we can solve $$\max( \sin(x)^2 + \sin(x^2) ), \quad \quad \lfloor x \rfloor = \mbox{prime}.$$

x = chebfun('x', [0 10]);
objective = chebfun( @(x) sin(x).^2 + sin(x.^2), [0 10] );
constrain = 0*x;
p = primes( 10 );
for j = 1:length(p)
constrain = constrain + ((p(j)<x) + (x<(p(j) + 1))-1);
end
g = objective.*constrain;
[mx, loc] = max( g );
plot(g,LW,lw), hold on, plot(loc, mx, 'r.', MS, ms), hold off
title(sprintf('constrained maximum = %1.3f',mx),FS,fs), set(gca,FS,fs)
ylim([-2 3]);

Here, to deal with the simple constraint on the independent variable a characteristic function was constructed and multiplied with the objective function.

### Another constraint

A similar thing can be done with other constraints. For example, for $x\in[0,10]$, consider

$$\min(\sin(x)^2+\sin(x^2)), \quad \quad |\sin(10x)|<{1\over 2}.$$

We can proceed like this:

x = chebfun('x', [0 10]);
objective = chebfun( @(x) sin(x).^2 + sin(x.^2), [0 10] );
constrain = (abs(sin(10*x)) < 1/2);
g = objective.*constrain;
[mx, loc] = max( g );
plot(g,LW,lw), hold on, plot(loc,mx,'r.',MS,ms), hold off
title(sprintf('constrained maximum = %1.3f', mx),FS,fs), set(gca,FS,fs)
ylim([-2 3]);

### Two-dimensional constrained optimization

Chebfun2 also has capabilities for global rootfinding and optimization and hence, also some constrained optimization. For example, consider maximizing

$$\cos((x-1/10)y)^2 + x \sin(3x+y)$$

over this heart-shaped region:

t = chebfun('t',[0,2*pi]);
x = 2*sin(t); y = 2*cos(t)-(1/2)*cos(2*t)-(1/4)*cos(3*t)-(1/8)*cos(4*t);
constrain = x + 1i*y;
plot(constrain,'k-',LW,lw), axis equal, title('Constraint',FS,fs)
axis([-3 3 -3 3])

Here is a contour plot of the objective function:

objective = chebfun2(@(x,y) cos((x-.1).*y).^2 + x.*sin(3*x+y), [-3 3 -3 3]);
contour(objective, LW, lw), axis equal, hold on
plot( constrain, 'k-', LW, lw), set(gca, FS, fs)
axis([-3 3 -3 3])

We can solve this by first finding all the local extrema of the objective function, restricting to those that lie inside the heart, and then picking the maximum:

r = roots( grad( objective ) );
idx = inpolygon(r(:,1), r(:,2), constrain );
r = r(idx,:);
plot(r(:,1), r(:,2), '.k', MS, ms)
max_inside = max(objective(r(:,1), r(:,2)));
max_boundary = max(objective(constrain));
max_overall = max(max_inside, max_boundary)

The maximum occurs inside the heart. Let's plot it:

[ignored, loc] = max(objective(r(:,1), r(:,2)));
plot( r(loc,1), r(loc,2), 'r.', MS, 40)
title(sprintf('Overall maximum = %1.3f', max_overall), FS, fs)

### Humble comment

The reason that Chebfun and Chebfun2 are able to perform some constrained optimization is by virtue of global optimization. For large scale problems this approach quickly becomes too computationally expensive, and there is a huge field of mathematics devoted to more efficient methods.

end

### Inpolygon

The script below is an overload of the inpolygon command in MATLAB. Functionality may be added into a new release of Chebfun, in which case this script will be removed.

function [in,on] = inpolygon(x,y,p)
%INPOLYGON True for points inside or on a region with chebfun boundary.
%   IN = INPOLYGON(X, Y, P) returns a matrix IN the size of X and Y. IN(j,k) = 1
%   if the point (X(j,k), Y(j,k)) is either strictly inside or on the edge of
%   the region with boundary given by the chebfun P; otherwise IN(j,k) = 0.
%
%   [IN ON] = INPOLYGON(X, Y, P) returns a second matrix, ON, which is the size
%   of X and Y. ON(j,k) = 1 if the point (X(j,k), Y(j,k)) is on the edge of the
%   polygonal region; otherwise ON(i,j) = 0.

tol = 100*eps;

p = p(:); % concatenate to column chebfun.
if ( all(abs(p(p.domain,'left') - p(p.domain,'right')) > 100*eps) )
error('CHEBFUN:inpolygon:closed','Chebfun must form closed curve.');
end

pgon = 1;   %Is it a polygon?
for k = 1:numel(p.funs{:})
if ( length(p.funs{k}) > 2 ), pgon = 0; break; end
end

% If it is a polygon use MATLAB inbuilt inpolygon command.
if ( pgon )
[in, on] = inpolygon( x, y, real(p.values{:}), imag(p.values{:}) );
return
end

% Not a polygon. Use Cauchy residue theorem (winding numbers).
in = zeros( numel(x), 1); on=in;
for k = 1:numel(x)
if (~isempty(roots(p - x(k) - 1i*y(k))))  % Does the curve & point intersect?
on(k) = 1;
end
end

left = 1:numel(x); left = left(~on);      % If not on, is it in?
% count winding number
for k = left
realp = real(p);
rts = roots(realp-x(k));              % Does curve cross vertical line?
if (isempty(rts)), continue, end      % No cross -> not in.
sgn = sign(imag(p(rts))-y(k));        % Are we above or below the point?
way = diff(realp);
way = sign(way(rts));                 % Clockwise or counterclockwise?
if ( abs(sum(sgn.*way)) > 1-tol ), in(k) = 1; end % Count up signs.
end

on = on | on;  % convert to boolean
in = in | on;

end
max_overall =
1.659510987079417