Python and Scipy Optimization implementation

Question

I wrote a code and I need some help about implementing an optmization method, maybe with scipy. If you can note below I have a variable "pD" that I need to vary in order to find a result for "abs(mFmin[i][j] - mReg[i][j]) > 1". mFmin, mReg and all the other calculations inside this while depends on "pD"

I wrote a poor code, just for testing, in order to vary "pD" comparing mFmin and mReg, but it's too slow and doesn't matter if I raise the error or not, that little code sucks.

I'm searching some optmization code in scipy library, but I can't see where I can implement this with my code. I think it's simples solving, but I have nobody to exchange ideas...

Note: pD is a matrix

Below, I attached the main part of the code:

for i in range(0,int(x)):
    for j in range(0,n):
        curso[i] = i*passo
        curso[0] = 0
        pD[i][j] = pZref
        mFmin[i][j] = 0
        mReg[i][j] = gref

# my doubt starts here
        while abs(mFmin[i][j] - mReg[i][j]) > 1:
            if mFmin[i][j] < mReg[i][j]:
                pD[i][j] = pD[i][j] + 0.0001
            else:
                pD[i][j] = pD[i][j] - 0.0001

            pZaux[i][j] = pE_*sqrt((pow(pZref/pE_,2)-pow(pA/pE_,2))*pow(mFmin[i][j]/gref,2)+pow(pA/pE_,2))

            vD[i][j] = pE_*vE_/pD[i][j]
            if pD[i][j]/pE_ > RPcr:
                psiR[i][j] = sqrt(pow(pD[i][j]/pE_,2/kappa)-pow(pD[i][j]/pE_,(kappa+1)/kappa))
            else:
                psiR[i][j] = psicrR
            if pZaux[i][j]/pD[i][j] > RPcr:
                psiF[i][j] = sqrt((2*kappa/(kappa-1))*(pow(pZaux[i][j]/pD[i][j],2/kappa)-pow(pZaux[i][j]/pD[i][j],(kappa+1)/kappa)))
            else:
                psiF[i][j] = psicrF

            mFmin[i][j] = 3600*psiF[i][j]*kmin*(fmin[j]/1000000)*sqrt(pD[i][j]*100000/vD[i][j])
            mReg[i][j] = 3600*psiR[i][j]*alpha*(fV[i][j]/1000000)*sqrt(2*kappa/(kappa-1)*(pE_*100000/vE_))

Thanks for reading!

MRM

First of all, use `numpy` and get rid of those pure python loops. — Andrey Sobolev, Mar 04 '15 at 18:38
Any tip? Actually, I'm using numpy so I can do the for's and all the math. I'm not a programmer experienced, I just know basic concepts...sorry... — Marcus R Marincek, Mar 04 '15 at 19:03

Balint Domokos · Answer 1 · 2015-03-04T20:54:29.420

You could take a look at the optimize package: http://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html

As a simple example let's say you are looking for the minimum of (x-3)**2. You define a function, which gets the input, and returns the function value. You pass this function to minimize, along with an initial guess x0.

import numpy as np
from scipy.optimize import minimize
def fn(x):
    return (3-x)**2
x0 = 0
res = minimize(fn, x0, method='nelder-mead', options={ 'disp': True})

This returns 3.0 as expected. You can define fn to have an input vector x, and then specify x0 as a vector the initial point with same dimension as expected by fn.

In the example 'nelder-mead' method is a simple algorithm which might have a long running time. If you know the grandient or your function to be minimized, you can use more advanced algorithms e.g. BFGS, and pass the gradient function as well, as described in the doc.

Thanks for the comment. I have already look the scipy docs, but I can't see how to put a method of optmization inside my program. — Marcus R Marincek, Mar 04 '15 at 19:05
I have updated my answer to provide an example and some explanation. — Balint Domokos, Mar 04 '15 at 20:56

score 0 · Answer 2 · answered Mar 05 '15 at 12:48

0

Andrey Sobolev suggests to use numpy because your if pD was a numpy array you would access its elements as pD[i,j] not pD[i][j], which is faster and simpler.

If I got it right, each [i,j] optimisation is independent of the other values of the array, right? Then, you just have to perform i*j optimisations (if the results were coupled, things would get more complicated...)

This problem may be too slow for 2 reasons. Because the optimizations slow, or because i and j are very large. In the later case you can probably accelerate it a lot using Numba, as the for loops are slow in Python.

As already suggested, use some method in scipy.optimize to do each of the optimizations. Be careful with the initial guess, specially if the function has more than one minimum.

answered Mar 05 '15 at 12:48

Ramon Crehuet

3,679
1
22
37

Andrey Sobolev, Balint Domokos and user2783943, thank you for your comments and interest. I will follow your instructions. As soon as I get it I put it here! Thank you all again! – Marcus R Marincek Mar 05 '15 at 21:21
user2783943, you are right, pD is independent. This program is about a valve lift where i is a lift value (it varies 0 mm to 20 mm for example) and j is the number of valves (it can calculate valve 1, valve 2, valve j...) as it varies we have an area fV and it gives a flow by an iteratively found pressure (pD). So as you can probably guess, each pD calculates one flow for each lift of each valve. – Marcus R Marincek Mar 05 '15 at 21:25

score 0 · Answer 3 · edited May 23 '17 at 10:30

0

I'd recommend importing numpy with:

import numpy as np

because your error it seems as your sqrt comes from the math module, and not the numpy module. This could be the cause of the error. See TypeError: only length-1 arrays can be converted to Python scalars while trying to exponentially fit data

Some more recommendations:

There is no need to specify 0 in slice ranges: psiR[0:,0:]==psiR[:,:]==psiR
create x0 with x0 = np.ones(n) * pZref if pZref is a scalar.
You should read some numpy manual to get familiar with these constructions.

edited May 23 '17 at 10:30

Community

1
1

answered Mar 07 '15 at 21:22

Ramon Crehuet

3,679
1
22
37

Ramon Crehuet, thanks for the advices! I will follow what you said. Have a nice week! – Marcus R Marincek Mar 08 '15 at 15:05

Python and Scipy Optimization implementation

3 Answers3