Using VBA to iteratively solve a system of equations

Question

I'm trying to solve the following system of equations in VBA:

I have looked for other similar questions and I can't find any with suitable solutions. I have already solved the equations in a worksheet by simply entering the equations as formulae is separate cells (which will initially create a circular reference warning) and the enabling iterative calculation) - given that I know DR will always be greater than PR, I thought an iterative approach of incrementally reducing PR from an initial value of DR would work in VBA. For reference, the worksheet used is shown below:

The formula in G9 is =(G8*B6)/(G10+273.15), i.e. equation 2 (see above), and the formula in G10 is =(B6+(B6*(-1+((G9*(B3-B35))/(B3-B35))^0.263)/B34))-B33*(B6+(B6*(-1+((G9*(B3-B35))/(B3-B35))^0.263)/B34)), i.e. equation 1 (see above).

When I try to do this programmatically in VBA by reducing PR incrementally from an initial value of DR, it doesn't work. My code is below:

Sub ChargeTempAndPressureCalculations()

Dim AP_hPa As Double
Dim AP_psi As Double
Dim TIn_C As Double
Dim TIn_K As Double
Dim PR As Double
Dim Ei As Double
Dim Et As Double
Dim Vci As Double
Dim DR As Double
Dim TOut_C As Double

AP_hPa = 1029 'Input
AP_psi = AP_hPa * 100 * 0.000145038

TIn_C = 15 'Input
TIn_K = TIn_C + 273.15

Et = 0.75 'Input
Ei = 0.75 'Input
Vci = 0.5 'Input
DR = 2.7103502887329 'Input
PR = DR

Do Until TOut_C = (TIn_K + (TIn_K * (-1 + ((PR * (AP_psi - Vci)) / (AP_psi - Vci)) ^ 0.263) / Et)) - Ei * (TIn_K + (TIn_K * (-1 + ((PR * (AP_psi - Vci)) / (AP_psi - Vci)) ^ 0.263) / Et)) And PR = (DR * TIn_K) / (TOut_C + 273)
PR = PR - 0.00000000001
Loop

Debug.Print "Charge air temperature = " & TOut_C
Debug.Print "Pressure Ratio = " & PR


End Sub

It's clearly the loop that's the issue but what is it that I'm doing wrong?

Edit:

I've split equation 1 to prevent the 'Equation too complex' error observed by another user, mentioned in the comments. I've also added a control to prevent the number of steps going above 1000.

Sub ChargeTempAndPressureCalculations()

Dim AP_hPa As Double
Dim AP_psi As Double
Dim TIn_C As Double
Dim TIn_K As Double
Dim PR As Double
Dim Ei As Double
Dim Et As Double
Dim Vci As Double
Dim DR As Double
Dim TOut_C As Double
Dim A As Double
Dim B As Double
Dim i As Integer

AP_hPa = 1029 'Input
AP_psi = AP_hPa * 100 * 0.000145038

TIn_C = 15 'Input
TIn_K = TIn_C + 273.15

Et = 0.75 'Input
Ei = 0.75 'Input
Vci = 0.5 'Input
DR = 2.7103502887329 'Input
PR = DR

Do Until i > 1000 Or (TOut_C = A - Ei * B And PR = (DR * TIn_K) / (TOut_C + 273))
'Spliting equation for TOut_C to simplify the expression and prevent an error
A = (TIn_K + (TIn_K * (-1 + ((PR * (AP_psi - Vci)) / (AP_psi - Vci)) ^ 0.263) / Et))
B = (TIn_K + (TIn_K * (-1 + ((PR * (AP_psi - Vci)) / (AP_psi - Vci)) ^ 0.263) / Et))

PR = PR - 0.00000000001
i = i + 1
Loop

Debug.Print "Charge air temperature = " & TOut_C
Debug.Print "Pressure Ratio = " & PR


End Sub

Having read the answer provided, I'm still none the wiser as to how to resolve my issue.

score 0 · Answer 1 · answered Mar 16 '17 at 10:01

0

Note that in:

Do Until TOut_C = ...

the = sign is a comparisson and not an assignment. As TOut_C has not been used yet, it is set to zero by VB and so you are comparing whether the right hand side is zero. This does not seem to be your intention, as you use TOut_C in the AND part as TOut_C + 273 which then would always be 273.

But if you want to compare with zero, then note that the RHS will probably never become zero in floating point arithmetic and you must compare with an "epsilon", a small value that is your precission threshold. For example:

Private Const eps = 0.00000000001   ' must be smaller than your step size

Do Until Abs(TOut_C - RHS) < eps

I leave fixing this to you. (I also get an error "Expression too complex" on my Excel version.)

answered Mar 16 '17 at 10:01

Paul Ogilvie

25,048
4
23
41

As you correctly surmised, it is not my intention to compare with zero. Given that's the case, I'm afraid your answer doesn't actually solve my problem. Is there more information I can provide to help you answer? – newishuser Mar 16 '17 at 10:54
You keep comparing the terms with `TOut_C` and `TOut_C` is zero. You now have moved some terms outside the loop, but they contain `PR` but `PR` is only updated in the loop. The terms will not be recomputed if they are outside the loop. – Paul Ogilvie Mar 16 '17 at 13:33
Could you post the excel sheet somewhere on the web so I can study it and see how you do the computation in the sheet? – Paul Ogilvie Mar 16 '17 at 13:36
I'm afraid I can't access any online storage services at the moment. I've expanded my explanation in the original question. In short, I just put the equations in separate cells, referencing each other (obviously this created a circular reference warning). I then enabled iterative calculation which solved the equations. Given that I know PR is always less than DR, reducing PR incrementally from a value of DR should eventually provide a valid solution. – newishuser Mar 16 '17 at 15:07
You seem to have a problem with understanding programming in general. I suggest you contact a colleague to help you. Without the sheet I can't help you. – Paul Ogilvie Mar 16 '17 at 15:36
You may well be right. I doubt my colleagues will be able to help but thank you for the suggestion. If you aren't completely fed up with my question yet, here is the [link](https://drive.google.com/open?id=0B_XwdXaxotUzT2Y2YVNKWVhYMXM) – newishuser Mar 16 '17 at 16:31
I've downloaded it; I will be back tomorrow. – Paul Ogilvie Mar 16 '17 at 17:33

score 0 · Answer 2 · edited May 23 '17 at 10:29

There is primarily two problems... In your Do..Loop you want the code to break on TOut_C = <something based on PR> AND PR = <something based on TOut_C>

Both equation however lead to a double datatype, this is almost imposibble to have that be an = comparison as the likelyhood of hitting that equal point are virtualy zero (as I explained here)

So, you would want to set a more flexible parameter such as TOut_C > 63 for example.

The second thing that I can find is that the equations are just that, equations. So they'll produce a number but as that number is the result of the equation, but what would define the endpoint ? TOut_C is never defined and subsequently compared as indicated, your code starts with PR being 2.71... but Tout_C is 0.

So could you elaborate more on the correlation between TOut_C and PR and what would be the solution values for either or both of them ? Or are you trying to solve a balance point i.e. for a which TOut_C do both functions achieve the same result ? (That would require mathematically re-writing either function to express the same output) So rewrite the function for PR to produce Tout_C...

The loop works (i.e. it loops) if you write it like so... It doesn't solve it, but at least you know it loops...

Do Until TOut_C = 15 And PR = 22
TOut_C = (TIn_K + (TIn_K * (-1 + ((PR * (AP_psi - Vci)) / (AP_psi - Vci)) ^ 0.263) / Et)) - Ei * (TIn_K + (TIn_K * (-1 + ((PR * (AP_psi - Vci)) / (AP_psi - Vci)) ^ 0.263) / Et))
PR = (DR * TIn_K) / (TOut_C + 273)
PR = PR - 0.00000000001
Loop

mtholen · Accepted Answer · 2017-03-17T01:35:30.900

I've gone over your calculation a bit more and I presume you want to stop if the increment of PR and TOut_C is very marginal...

The below code does exactly that. It calculates TOut_C at a given PR, it subsequently calculates the PR corresponding to that TOut_C, the loop then substitutes the newly calculated PR into the TOut_C calculation and so on.

It calculates the difference between the substitute calculation and the previous calculation and if there is no longer a 'large' offset between the two it stops the loop.

Sub ChargeTempAndPressureCalculations()

Dim AP_hPa As Double
Dim AP_psi As Double
Dim TIn_C As Double
Dim TIn_K As Double
Dim PR As Double
Dim Ei As Double
Dim Et As Double
Dim Vci As Double
Dim DR As Double
Dim TOut_C As Double

AP_hPa = 1029 'Input
AP_psi = AP_hPa * 100 * 0.000145038

TIn_C = 15 'Input
TIn_K = TIn_C + 273.15

Et = 0.75 'Input
Ei = 0.75 'Input
Vci = 0.5 'Input
DR = 2.7103502887329 'Input
PR = DR

dTOut_C = 1     'Set to arbitrary number to initialize the loop
dPR = 1         'Set to arbitrary number to initialize the loop

Do Until dPR < 0.0000000001 And dTOut_C < 0.0000000001
'Calculate the TOut_C and PR
TOut_C = (TIn_K + (TIn_K * (-1 + ((PR * (AP_psi - Vci)) / (AP_psi - Vci)) ^ 0.263) / Et)) - Ei * (TIn_K + (TIn_K * (-1 + ((PR * (AP_psi - Vci)) / (AP_psi - Vci)) ^ 0.263) / Et))
PR = (DR * TIn_K) / (TOut_C + 273)

'Calculate difference relative to last calculation
dPR = PR - PR0
dTOut_C = TOut_C - TOut_C0

'Set the last calculation as previous calculation and re-do loop
PR0 = PR
TOut_C0 = TOut_C

Loop

Debug.Print "Charge air temperature = " & TOut_C
Debug.Print "Pressure Ratio = " & PR

End Sub

Given your initial input the output is:

Charge air temperature = 93.2076926574912
Pressure Ratio = 2.13263525413933

Was that what you were looking for ?

PS: Technically what you should do is mathematically rewrite the equation to solve for PR based on DR as that is in essence what you are doing...

Many thanks. I realise that rearranging the equations would have been a more elegant solution but using iterative calculation in the worksheet was a shortcut and it got me wondering whether a similar approach could be used in VBA. You've satisfied my curiosity. Thanks again. — newishuser, Mar 17 '17 at 07:53

score 0 · Answer 4 · answered Mar 17 '17 at 09:27

See http://www.decisionmodels.com/calcsecretsc.htm on how Excel processes circular references. Basically, it just calculates each cell, ignoring the circular reference and then it updates the values in each iteration.

Applying this to your VBA routine produces the following sub routine:

Sub ChargeTempAndPressureCalculations()

    ' input variables
    Dim AP_hPa As Double
    Dim AP_psi As Double
    Dim TIn_C As Double
    Dim TIn_K As Double
    Dim Ei As Double
    Dim Et As Double
    Dim Vci As Double
    Dim DR As Double

    ' temporary variables
    Dim Td As Double
    Dim Pd As Double
    Dim A As Double

    ' output variables
    Dim TOut_C As Double
    Dim PR As Double

    ' iteration control
    Dim eps As Double
    Dim i As Integer

    AP_hPa = 1029
    AP_psi = AP_hPa * 100 * 0.000145038

    TIn_C = 15
    TIn_K = TIn_C + 273.15

    Et = 0.75
    Ei = 0.75
    Vci = 0.5
    DR = 2.7103502887329
    PR = DR

    eps = 0.00000000001
    i = 0
    Do
        Td = TOut_C     ' remember values from previous iteration ( 'd' means 'delta')
        Pd = PR

        A = (TIn_K + (TIn_K * (-1 + ((PR * (AP_psi - Vci)) / (AP_psi - Vci)) ^ 0.263) / Et))
        TOut_C = A - Ei * A
        PR = (DR * TIn_K) / (TOut_C + 273)
        i = i + 1

        Debug.Print TOut_C & ", " & PR & "(" & Abs(Td - TOut_C) & ", " & Abs(Pd - PR) & ")" ' show progression

        ' loop until the difference is less than eps or max iterations reached
    Loop While (i < 100 And (Abs(Td - TOut_C) > eps And Abs(Pd - PR) > eps))


    Debug.Print "Charge air temperature = " & TOut_C
    Debug.Print "Pressure Ratio = " & PR
    Debug.Print "number of iterations: " & i

End Sub

Output:

100.835921416446, 2.08911822261291(100.835921416446, 0.621232066119993)
92.5738330344343, 2.13633297880164(8.26208838201211, 4.72147561887288E-02)
93.2611120407512, 2.13232420812257(0.687279006316899, 4.00877067907057E-03)
93.2031960023579, 2.13266144103602(5.79160383933299E-02, 3.37232913455665E-04)
93.2080712078647, 2.13263304962791(4.87520550686327E-03, 2.83914081067316E-05)
93.2076607895546, 2.13263543972443(4.10418310181626E-04, 2.39009651137323E-06)
93.2076953402811, 2.13263523851592(3.45507265677725E-05, 2.01208508077144E-07)
93.2076924316548, 2.1326352554545(2.90862628560262E-06, 1.6938581648418E-08)
93.2076926765153, 2.13263525402854(2.44860444809092E-07, 1.42596112695514E-09)
93.2076926559019, 2.13263525414858(2.06133563551703E-08, 1.200430865822E-10)
93.2076926576372, 2.13263525413848(1.73530168012803E-09, 1.01052499701382E-11)
93.2076926574912, 2.13263525413933(1.46073375617561E-10, 8.5043083686287E-13)
Charge air temperature = 93.2076926574912
Pressure Ratio = 2.13263525413933
number of iterations: 12

Using VBA to iteratively solve a system of equations

Edit:

4 Answers4