Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Optimization using Heaviside function

Subject: Optimization using Heaviside function

From: Petros

Date: 13 Jun, 2010 17:37:03

Message: 1 of 18

I try to run a one-variable optimization of a function that makes use of the heaviside function.

I have:
Q={700 if x(1)<5, 600 if 5<x(1)<10, 200 if 10<x(1)} or using heaviside functions,
Q=700-100*heaviside(x(1)-5)-400*heaviside(x(1)-10)

also I have
P=100+100*heaviside(x(1)-6)+400*heaviside(x(1)-13)

bound: x(1)>0

I want to minimize f=Q-P

I always get:
Optimization completed because the objective function is non-decreasing in
feasible directions, to within the default value of the function tolerance,
and constraints were satisfied to within the default value of the constraint tolerance.

and I get as solution the point I give as starting point. I believe the problem is the heavisides?

Thanks.

Subject: Optimization using Heaviside function

From: Roger Stafford

Date: 13 Jun, 2010 18:25:06

Message: 2 of 18

"Petros " <p3tris@gmail.com> wrote in message <hv34vv$91u$1@fred.mathworks.com>...
> I try to run a one-variable optimization of a function that makes use of the heaviside function.
>
> I have:
> Q={700 if x(1)<5, 600 if 5<x(1)<10, 200 if 10<x(1)} or using heaviside functions,
> Q=700-100*heaviside(x(1)-5)-400*heaviside(x(1)-10)
>
> also I have
> P=100+100*heaviside(x(1)-6)+400*heaviside(x(1)-13)
>
> bound: x(1)>0
>
> I want to minimize f=Q-P
>
> I always get:
> Optimization completed because the objective function is non-decreasing in
> feasible directions, to within the default value of the function tolerance,
> and constraints were satisfied to within the default value of the constraint tolerance.
>
> and I get as solution the point I give as starting point. I believe the problem is the heavisides?
>
> Thanks.

  It is not appropriate to use the standard optimization routines for step functions like this. These optimization routines typically depend heavily on their input functions being continuous and preferably possessing continuous derivatives, which is most certainly not true of the heavyside function.

  In any case you are making a very hard job out of an easy one. Just evaluate P-Q at and on either side of each step and find the minimum of these. What could be simpler?

Roger Stafford

Subject: Optimization using Heaviside function

From: Walter Roberson

Date: 13 Jun, 2010 18:33:27

Message: 3 of 18

Petros wrote:
> I try to run a one-variable optimization of a function that makes use of
> the heaviside function.

> I always get:
> Optimization completed because the objective function is non-decreasing
> in feasible directions, to within the default value of the function
> tolerance,
> and constraints were satisfied to within the default value of the
> constraint tolerance.

Which minimizer are you using? A number of the minimizers assume that
the differential of the function is continuous, which is not the case
for Heaviside. Or to put it another way, your graph is too smooth
(indeed, it is mostly flat), and most minimizers need a slope to work with.

You don't need a full minimizer for functions such as that: you only
have to test both sides of each boundary condition.

Subject: Optimization using Heaviside function

From: Petros

Date: 13 Jun, 2010 18:43:03

Message: 4 of 18

Thanks for replying. Really the function is not that one. The real function is a 2-variable function with quadratic parts and heaviside functions. I just wanted to start with something simpler. I cannot avoid the heaviside functions unfortunately, because they represent decisions made by a controller. So, the functions are

P={a if 0<x1<5; b if 5<x1} and
Q={c if 0<x2<19; d if 19<x2} and
L=(Q-P-Q0)*FINE(Q-P-Q0) where
FINE(Q-P-Q0)={e if 0<(Q-P-Q0)<100; f if 100<(Q-P-Q0)<200}

And I try to maximize (Q-P-L) or minimize -(Q-P-L).
It's not something I can do by hand for repeated problems and I thought matlab could help me out.

Thanks again.

Subject: Optimization using Heaviside function

From: Petros

Date: 13 Jun, 2010 19:06:03

Message: 5 of 18

If I use an approximation like: heaviside(x)=1/(1+exp(-50000*x)) which gives a good approximation of the heaviside funtion, will it be ok for using fmincon?

Subject: Optimization using Heaviside function

From: Roger Stafford

Date: 13 Jun, 2010 19:48:03

Message: 6 of 18

"Petros " <p3tris@gmail.com> wrote in message <hv38rn$sfv$1@fred.mathworks.com>...
> Thanks for replying. Really the function is not that one. The real function is a 2-variable function with quadratic parts and heaviside functions. I just wanted to start with something simpler. I cannot avoid the heaviside functions unfortunately, because they represent decisions made by a controller. So, the functions are
>
> P={a if 0<x1<5; b if 5<x1} and
> Q={c if 0<x2<19; d if 19<x2} and
> L=(Q-P-Q0)*FINE(Q-P-Q0) where
> FINE(Q-P-Q0)={e if 0<(Q-P-Q0)<100; f if 100<(Q-P-Q0)<200}
>
> And I try to maximize (Q-P-L) or minimize -(Q-P-L).
> It's not something I can do by hand for repeated problems and I thought matlab could help me out.
>
> Thanks again.

  You haven't stated what Q0 is, but if it's like P and Q, then you still have a step function, though in two variables. That is, the function is constant over various rectangular regions in the x-y plane, (or I should say x1-x2 plane.) In spite of the '*' operation it is in no sense a quadratic function. The standard optimization routines still cannot handle these flat planar regions properly.

  I think you have no choice but to evaluate the function on either side of its discontinuities and select the minimum or maximum. After all, the answer lies among only finitely many possible sums, differences, and products of the numbers, a, b, c, d, e, and f. It is a mistake to invoke routines that are attuned to deal with a (nearly) infinite continuum of values?

Roger Stafford

Subject: Optimization using Heaviside function

From: Roger Stafford

Date: 13 Jun, 2010 20:25:04

Message: 7 of 18

  (Ignoring the as yet unspecified Q0,) you have discontinuities along the x1 axis only at 0 and 5, and along the x2 axis at 0 and 19. Define

 x1 = [0-1,(0+5)/2,5+1];
 x2 = [0-1,(0+19)/2,19+1];

Now evaluate Q-P-L at each of the nine combinations of values in x1 and x2. That is, at

 [X1,X2] = meshgrid(x1,x2);

and select the maximum of these nine Q-P-L values.

  As long as your function consists only of constant rectangular regions as in the case you have described, this kind of procedure should be all that is necessary to find the maximum. It is only when the total number of combined discontinuities are too large for practical computation that more heroic methods become necessary.

Roger Stafford

Subject: Optimization using Heaviside function

From: Petros

Date: 13 Jun, 2010 20:25:05

Message: 8 of 18

I understand what you are saying. But, these functions get really complicated. How do I check finite number of possible solutions? I have managed to plot a few of them with syms and ezsurf() after that what?

Q0 is a constant and the equation sometimes goes:

P={c*x2^2 if x2<5; d*x2^2+e*x2 if x2>5}

that's why I said quadratic. If I use the approximation I stated above and use the fmincon, will it be ok?

Subject: Optimization using Heaviside function

From: John D'Errico

Date: 13 Jun, 2010 20:59:04

Message: 9 of 18

"Petros " <p3tris@gmail.com> wrote in message <hv3a6r$ip3$1@fred.mathworks.com>...
> If I use an approximation like: heaviside(x)=1/(1+exp(-50000*x)) which gives a good approximation of the heaviside funtion, will it be ok for using fmincon?

NO!!!!!!!

You don't understand that the problem is the non-differentiability
of your function, when you then try to use an optimizer that is
designed to work on a differentiable objective function.

Creating a function that is IDENTICAL to a step function in ALL
respects when used in floating point arithmetic is not a valid
solution, since that function is also essentially non-differentiable
in terms of floating point arithmetic.

You can't fool mathematics by being tricky.

John

Subject: Optimization using Heaviside function

From: Roger Stafford

Date: 13 Jun, 2010 22:42:04

Message: 10 of 18

"Petros " <p3tris@gmail.com> wrote in message <hv3er0$nk7$1@fred.mathworks.com>...
> I understand what you are saying. But, these functions get really complicated. How do I check finite number of possible solutions? I have managed to plot a few of them with syms and ezsurf() after that what?
>
> Q0 is a constant and the equation sometimes goes:
>
> P={c*x2^2 if x2<5; d*x2^2+e*x2 if x2>5}
>
> that's why I said quadratic. If I use the approximation I stated above and use the fmincon, will it be ok?

  I have bad news for you Petros! As John has pointed out, using exp(-50000*x) will do you no good at all. It is necessary to face the discontinuities as they stand and deal with them effectively.

  It would be easily possible to use combinations of heavyside-type equations in such a way that would produce many not necessarily rectangular regions in the x1-x2 space in which the boundaries are discontinuous for the objective function and inside of which that function is variable in complex ways even though "smooth". If you contemplate something of this kind, I believe there is no option for you but to handle each such region separately, using an optimization routine on each if necessary. I know of no matlab routine that would automatically carry out such an analysis for you.

  Wasn't it Euclid who said to King Ptolemy, "there is no Royal Road to geometry"? An analogous remark might conceivably be made of functions composed of complicated heavyside-type combinations.

Roger Stafford

Subject: Optimization using Heaviside function

From: Greg Heath

Date: 14 Jun, 2010 02:06:44

Message: 11 of 18

On Jun 13, 1:37 pm, "Petros " <p3t...@gmail.com> wrote:
> I try to run a one-variable optimization of a function that makes use of the heaviside function.
>
> I have:
> Q={700 if x(1)<5, 600 if 5<x(1)<10, 200 if 10<x(1)} or using heaviside functions,
> Q=700-100*heaviside(x(1)-5)-400*heaviside(x(1)-10)
>
> also I have
> P=100+100*heaviside(x(1)-6)+400*heaviside(x(1)-13)
>
> bound: x(1)>0
>
> I want to minimize f=Q-P
>
> I always get:
> Optimization completed because the objective function is non-decreasing in
> feasible directions, to within the default value of the function tolerance,
> and constraints were satisfied to within the default value of the constraint tolerance.
>
> and I get as solution the point I give as starting point. I believe the problem is the heavisides?
>
> Thanks.

If that, indeed, is your ONLY problem then you probably can
sneak up on the answer by replacing the Heaviside by a
sequence of logsigs (if you have the NN Toolbox otherwise
just use the definition)

logsig(a(n)*x) = 1/(1+exp(-a(n)*x))

a(n) > a(n-1) >= a(1) = 1

Hope this helps.

Greg

Subject: Optimization using Heaviside function

From: Petros

Date: 14 Jun, 2010 07:58:03

Message: 12 of 18

Greg Heath <heath@alumni.brown.edu> wrote in message <dde7f674-e072-434d-8213-10ee9a1a9d2e@5g2000vbf.googlegroups.com>...
> On Jun 13, 1:37 pm, "Petros " <p3t...@gmail.com> wrote:
> > I try to run a one-variable optimization of a function that makes use of the heaviside function.
> >
> > I have:
> > Q={700 if x(1)<5, 600 if 5<x(1)<10, 200 if 10<x(1)} or using heaviside functions,
> > Q=700-100*heaviside(x(1)-5)-400*heaviside(x(1)-10)
> >
> > also I have
> > P=100+100*heaviside(x(1)-6)+400*heaviside(x(1)-13)
> >
> > bound: x(1)>0
> >
> > I want to minimize f=Q-P
> >
> > I always get:
> > Optimization completed because the objective function is non-decreasing in
> > feasible directions, to within the default value of the function tolerance,
> > and constraints were satisfied to within the default value of the constraint tolerance.
> >
> > and I get as solution the point I give as starting point. I believe the problem is the heavisides?
> >
> > Thanks.
>
> If that, indeed, is your ONLY problem then you probably can
> sneak up on the answer by replacing the Heaviside by a
> sequence of logsigs (if you have the NN Toolbox otherwise
> just use the definition)
>
> logsig(a(n)*x) = 1/(1+exp(-a(n)*x))
>
> a(n) > a(n-1) >= a(1) = 1
>
> Hope this helps.
>
> Greg

Sorry, didn't get that NN thing.

I understand what all you people are saying and I agree. The problem is that with a rough estimation I have made the full scale problem will have more than 100 000 "pieces". That is why I hoped for the optimization. I can do it brute force (with x1 x2 increments of 0.1) but... Normally, if I aggregate all the heaviside in large number I can approximate each function with fit curve 2a+b*x+c*x^2+..." and go from there. I will go with the brute force program now and I will see how it goes.

Thanks for help. You gave me plenty ideas and helped me focus.

Subject: Optimization using Heaviside function

From: Greg Heath

Date: 14 Jun, 2010 13:10:57

Message: 13 of 18

On Jun 14, 3:58 am, "Petros " <p3t...@gmail.com> wrote:
> Greg Heath <he...@alumni.brown.edu> wrote in message <dde7f674-e072-434d-8213-10ee9a1a9...@5g2000vbf.googlegroups.com>...
> > On Jun 13, 1:37 pm, "Petros " <p3t...@gmail.com> wrote:
> > > I try to run a one-variable optimization of a function that makes use of the heaviside function.
>
> > > I have:
> > > Q={700 if x(1)<5, 600 if 5<x(1)<10, 200 if 10<x(1)} or using heaviside functions,
> > > Q=700-100*heaviside(x(1)-5)-400*heaviside(x(1)-10)
>
> > > also I have
> > > P=100+100*heaviside(x(1)-6)+400*heaviside(x(1)-13)
>
> > > bound: x(1)>0
>
> > > I want to minimize f=Q-P
>
> > > I always get:
> > > Optimization completed because the objective function is non-decreasing in
> > > feasible directions, to within the default value of the function tolerance,
> > > and constraints were satisfied to within the default value of the constraint tolerance.
>
> > > and I get as solution the point I give as starting point. I believe the problem is the heavisides?
>
> > > Thanks.
>
> > If that, indeed, is your ONLY problem then you probably can
> > sneak up on the answer by replacing the Heaviside by a
> > sequence of logsigs (if you have the NN Toolbox otherwise
> > just use the definition)
>
> > logsig(a(n)*x) = 1/(1+exp(-a(n)*x))
>
> > a(n) > a(n-1) >= a(1) = 1
>
> > Hope this  helps.
>
> > Greg
>
> Sorry, didn't get that NN thing.

If you have the NN Tbx, then you can use the function "logsig".
Otherwise, just use the definition and "exp".


Greg

Subject: Optimization using Heaviside function

From: Petros

Date: 14 Jun, 2010 13:24:04

Message: 14 of 18

I do have the NN toolbox (because I work at my university's pc which has all toolbox) but even if I use that, how do I solve the problem later?

Subject: Optimization using Heaviside function

From: Steven Lord

Date: 14 Jun, 2010 14:20:34

Message: 15 of 18


"Petros " <p3tris@gmail.com> wrote in message
news:hv34vv$91u$1@fred.mathworks.com...
>I try to run a one-variable optimization of a function that makes use of
>the heaviside function.
>
> I have:
> Q={700 if x(1)<5, 600 if 5<x(1)<10, 200 if 10<x(1)} or using heaviside
> functions,
> Q=700-100*heaviside(x(1)-5)-400*heaviside(x(1)-10)
>
> also I have P=100+100*heaviside(x(1)-6)+400*heaviside(x(1)-13)
>
> bound: x(1)>0
>
> I want to minimize f=Q-P

So Q's behavior changes at x1 = 5 and x1 = 10, while P's behavior changes at
x1 = 6 and x1 = 13.

Sounds like you have several problems here:

f(1) = minimize Q-P on the interval (-Inf, 5)
f(2) = minimize Q-P on the interval (5, 6)
f(3) = minimize Q-P on the interval (6, 10)
f(4) = minimize Q-P on the interval (10, 13)
f(5) = minimize Q-P on the interval (13, Inf)

Now the minimum of your whole function is the minimum value among the f's.

*snip*

> and I get as solution the point I give as starting point. I believe the
> problem is the heavisides?

The problem is the discontinuitites, as Roger and John have said.

--
Steve Lord
slord@mathworks.com
comp.soft-sys.matlab (CSSM) FAQ: http://matlabwiki.mathworks.com/MATLAB_FAQ
To contact Technical Support use the Contact Us link on
http://www.mathworks.com

Subject: Optimization using Heaviside function

From: Greg Heath

Date: 15 Jun, 2010 17:14:36

Message: 16 of 18

On Jun 14, 9:24 am, "Petros " <p3t...@gmail.com> wrote:
> I do have the NN toolbox (because I work at my university's pc which has all toolbox) but even if I use that, how do I solve the problem later?

I thought you said that the discontinuous heaviside was
your only problem.If that is true, try to construct a
sequence of solutions where each new solution is based
on replacing heaviside(z) with logsig(a(n)*z) with
a(n) > a(n-1) >= 1.

Hope this helps.

Greg

Subject: Optimization using Heaviside function

From: Johan Löfberg

Date: 16 Jun, 2010 10:33:08

Message: 17 of 18

Greg Heath <heath@alumni.brown.edu> wrote in message <85578c57-1844-440e-952c-6dda746ad5c3@t10g2000yqg.googlegroups.com>...
> On Jun 14, 9:24 am, "Petros " <p3t...@gmail.com> wrote:
> > I do have the NN toolbox (because I work at my university's pc which has all toolbox) but even if I use that, how do I solve the problem later?
>
> I thought you said that the discontinuous heaviside was
> your only problem.If that is true, try to construct a
> sequence of solutions where each new solution is based
> on replacing heaviside(z) with logsig(a(n)*z) with
> a(n) > a(n-1) >= 1.
>
> Hope this helps.
>
> Greg


If you are willing to use a mixed-integer solver (glpk, gurobi,cplex, etc), the following code implements my interpretation of your problem, in the modelling language YALMIP (free matlab toolbox)

d1 = sdpvar(2,1);
d2 = binvar(2,1);
d3 = binvar(2,1);

sdpvar L Q P x1 x2

a = 1;
b = 2;
c = 3;
d = 4;
e = 5;
f = 6;
Q0 = 7;

Constraints1 = [implies(d1(1),[P==a, 0<x1<5]),implies(d1(2),[P==b, x1>5]),sum(d1)==1]
Constraints2 = [implies(d2(1),[Q==c, 0<x2<19]),implies(d2(2),[Q==d, x2>19]),sum(d2)==1];
Constraints3 = [implies(d3(1),[L == (Q-P-Q0)*e, 0<Q-P-Q0<100])];
Constraints4 = [implies(d3(2),[L == (Q-P-Q0)*f, 100<Q-P-Q0<200]),sum(d3)==1]
Constraints = [Constraints1,Constraints2,Constraints3,Constraints4,-1000<[L Q P x1 x2]<1000];
solvesdp(Constraints,L)

Subject: Optimization using Heaviside function

From: Johan Löfberg

Date: 16 Jun, 2010 10:36:07

Message: 18 of 18

whoops, should be

d1 =binvar(2,1);



"Johan Löfberg" <loefberg@control.ee.ethz.ch> wrote in message <hva994$jg3$1@fred.mathworks.com>...
> Greg Heath <heath@alumni.brown.edu> wrote in message <85578c57-1844-440e-952c-6dda746ad5c3@t10g2000yqg.googlegroups.com>...
> > On Jun 14, 9:24 am, "Petros " <p3t...@gmail.com> wrote:
> > > I do have the NN toolbox (because I work at my university's pc which has all toolbox) but even if I use that, how do I solve the problem later?
> >
> > I thought you said that the discontinuous heaviside was
> > your only problem.If that is true, try to construct a
> > sequence of solutions where each new solution is based
> > on replacing heaviside(z) with logsig(a(n)*z) with
> > a(n) > a(n-1) >= 1.
> >
> > Hope this helps.
> >
> > Greg
>
>
> If you are willing to use a mixed-integer solver (glpk, gurobi,cplex, etc), the following code implements my interpretation of your problem, in the modelling language YALMIP (free matlab toolbox)
>
> d1 = sdpvar(2,1);
> d2 = binvar(2,1);
> d3 = binvar(2,1);
>
> sdpvar L Q P x1 x2
>
> a = 1;
> b = 2;
> c = 3;
> d = 4;
> e = 5;
> f = 6;
> Q0 = 7;
>
> Constraints1 = [implies(d1(1),[P==a, 0<x1<5]),implies(d1(2),[P==b, x1>5]),sum(d1)==1]
> Constraints2 = [implies(d2(1),[Q==c, 0<x2<19]),implies(d2(2),[Q==d, x2>19]),sum(d2)==1];
> Constraints3 = [implies(d3(1),[L == (Q-P-Q0)*e, 0<Q-P-Q0<100])];
> Constraints4 = [implies(d3(2),[L == (Q-P-Q0)*f, 100<Q-P-Q0<200]),sum(d3)==1]
> Constraints = [Constraints1,Constraints2,Constraints3,Constraints4,-1000<[L Q P x1 x2]<1000];
> solvesdp(Constraints,L)

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us