scipy.optimize with multiple bounds, constraints and continuous fields

Question

I want to optimize the operation of a CHP plant over a requested power profile. Therefore I define a power profile, which should be followed by the CHP plant as much as possible. Multiple bounds and constraints must be applied to represent the realistic operation of a CHP plant. This for example includes that the CHP can bei either on or off and while on, its power modulation can only be set to a specific percentage range.

Here is a minimum working example with short explanations:

import scipy.optimize as opt
import numpy as np

x = np.arange(200)  # dummy x vector
poly_profile = np.array(  # 7th degree polynome fit of profile
    [-2.14104340e-11,  1.85108903e-08, -6.66697810e-06,  1.29239710e-03,
     -1.45110876e-01,  9.40324129e+00, -3.24548750e+02,  4.60006330e+03])
poly_fun = np.poly1d(poly_profile)  # make poly fun
profile = poly_fun(x[65:196])
x0 = np.zeros_like(profile)  # all zeros as starting values

def optifun(x, profile):  # define minimization fun
    return - np.sum(profile * x)

bnds_hi = opt.Bounds(0.3, 1)  # upper bounds
bnds_lo = opt.Bounds(0, 0)  # lower bounds

res = opt.minimize(
    optifun, x0, args=(profile), bounds=bnds_hi,
    constraints={'type': 'eq', 'fun': lambda x:  np.sum(x*40) - 2000},
    method='SLSQP')
plt.plot(res.x)
plt.plot(profile)

So these are the bounds I want to use:

• (x == 0) or (0.3 <= x <= 1), for any value in the array x
  This means that x, the modulation degree of the total CHP power, can be either 0 (switched off) or >0.3 and <= 1. But I can either specify the lower bounds OR the upper bounds. Only specifying the upper bounds makes it impossible to "switch the CHP off", while setting the lower bounds to bnds_lo = opt.Bounds(0, 1)
  will enable the CHP plant to operate in a non-realistic operation point (between 0% and 30% of power modulation).
  Is there any way to make this work with the bounds as specified in the minimum working example? Specifically: Can I set two kinds of bounds at the same time, like bounds=[bnds_lo, bnds_hi]?
  I guess this is a mixed integer linear programming problem, but shouldn't COBYLA or SLSQP be able to handle this? If not: Is there any workaround?

And the constraints I want to use:

• np.sum(x*40) - 450
  Limits the thermal output to some thermal storage capacity. Here 40 is the thermal output power while 450 is the remaining storage capacity. This is fairly easy to implement.

• Limit the number of starts of the CHP plant. As an example let's assume

  bnds_lo = opt.Bounds(0, 1)  # lower bounds
  res = opt.minimize(
      optifun, x0, args=(profile), bounds=bnds_lo,
      constraints={'type': 'eq', 'fun': lambda x:  np.sum(x*40) - 1000},
      method='SLSQP')
  

  This leads to 3 periods of CHP plant operation. Is there any way to limit this? I was thinking about adding a specific constraints function which counts the positive diffs after a leading 0, but I was not able to make anything like this work (for example since most x are not exactly 0, since bounds are set to (0, 1). But other issues may also be the reason)...

• Set a minimum consecutive run time of the CHP plant. That means that it should be favourable to have at least, let's say 5, consecutive x != 0. I thought about trying something similar as in my last point (limiting the number of starts), but was also not able to work out something useful. This is by far the least important problem.

To solve these issues I also tried using scipy.optimize.LinearConstraings and NonlinearConstraings But method='trust-constr' requires a jac (as far as I read on github this seems to be a bug) and thus I wasn't able to make it work.

Is there any way I can make this work? Especially specifying multiple bounds is important.

Thanks in advance!

Sincerely, Scotty

Answer 1

profile * x0 in your code gives
"ValueError: operands could not be broadcast together with shapes (131,) (200,)".

Just guessing, is x_t a product onoff_t * xon_t
with onoff_t = 0 or 1
and 0.3 <= xon_t <= 1 at each t in 0 .. T ?
I.e. for T = 5 there are 2^5 possible onoff sequences, 00000 00001 00010 .. 11111 ?

If so, maximizng sum 0:T w_t * onoff_t * xon_t with a fixed weight function w_t is trivial:
where w_t <= 0: onoff_t = 0, off
where w_t > 0: onoff_t = 1, on, and xon_t = 1, max.
So that can't be your question -- please clarify.

If onoff_t is further constrained to switch only twice, 0... 1... 0..., then the number of possible sequences is small enough to just try them all, along the lines:

def pulse_generator( T=200, minwidth=5 ):
    """ -> arrays of T floats, 0... 1... 0... """
    for t0 in xrange( 1, T ):
        for t1 in xrange( t0 + minwidth, T ):
            pulse = np.zeros( T )
            pulse[t0:t1] = 1
            yield pulse

for pulse in pulse_generator( T ):
    print "pulse:", pulse
    optimize myfunction( pulse * xon ), 0.3 <= xon <= 1

Switching 4 times, 0... 1... 0... 1... 0..., is similar. (How many such pulses are there for a given T ? See wikipedia Stars and bars -- amazing.)

---

Added: I'm no expert, but isn't on-off aka bang-bang control very sensitive to tiny changes, a bit earlier or later ? A program(mer) can spend a lot of time dithering, down in the noise. How about 2 phases, coarse-grid then fine-grid --

1. split the time 0:T into say 10 pieces, run all 2^10 = 1024 on-off sequences
   1a. look at the best ones closely -- any pattern ?
2. move their edges by half steps, T / 20.

See also: google "discrete optimization" multigrid ... and Grid search.