Matlab fsolve


Hey everyone, long time no post :). It’s also been some time since I have been working in Matlab, but today I found enough time to waste. I am working on a UNIQUAC model for predicting binary solution activities, and also have been trying to generate a X-Y plot for distillation. I ran into a problem however, I have a function:

%Uniquac.m
function gamma = Uniquac(r_matrix, q_matrix, a_matrix, T, x1)
... (More code here)

The extra parameters in the function above are really nice, I might have to call several different binary components with the same function. While I could hard-code the constants r_matrix, q_matrix, and a_matrix… I really feel cheap doing so. Why does it matter so much? I don’t want to corrode my perfectly general code with very specific constants. Consider the very real possibility that you might want to use this function in fsolve, a numeric solver in Matlab. Fsolve requires:

1. That the first input be a function.
2. The second input must be a ‘guess’ for a possible solution.

The problem is that when you use a function like Uniquac.m, fsolve(@Uniquac, x0) says to Matlab: “Use the matrix x0 as the starting point for meeting my constraints.” But what is x0?

Well, it’s the inputs of Uniquac.m… which are r_matrix, q_matrix, a_matrix, T and x1. But hold on… maybe I only want to allow T to change and if r_matrix, q_matrix, a_matrix, and x1 are supposed to be constants, then it would be more than desirable if Matlab didn’t change them. But in this configuration Matlab cannot distinguish between the ‘constants’ and the variables that you want to change.

So what then? How to fix this problem? Well, you could go back and ruin all that beautiful code you just wrote by hard-coding constants and making global variables… but you could also go back to dragging your fists and beating things with a club.

OR we could do something VERY hip and gets rid of the bad vibes. For this we pull out a seemingly useless friend: the anonymous function. We covered these previously but here’s the angle:

%Some other script than Uniquac.m
shortUniquac = @(T) Uniquac(r_matrix, q_matrix, a_matrix, T, x1)

Yea! Now you’ve got a good as new function shortUniquac(T) that depends only on T… without a bunch of band-aids. Now for fsolve:

%fsolver
answer = fsolve(shortUniquac, T0)

If you would like the files for the Uniquac stuff, leave a comment with your email (do not put “example@gmail.com”, rather “example \ gmail” to avoid spam)  and what you intend to use it for and I’ll send you a link.

——-Updated 04/23/2012———

User List

  • Changsoo Lee — Yonsei Univ. Korea
  • Bobby J  — The University of Edinburgh
  • Santiago — University of Los Andes; Bogotá, Colombia
Advertisements

About pavlovsr

I'm a college student, just here to get some stuff out. Hope this helps. View all posts by pavlovsr

17 responses to “Matlab fsolve

  • aidin

    I would like mfile for unifac model

    • pavlovsr

      Aidin,

      Because I release the Unifac model as a non-profit, non-commercial use application I require that you send me an email that is associated with your learning institution (preferably a .edu) and what your intended use for the program is. Provided with these I can release the code to you, please keep in mind that it comes with ABSOLUTELY NO warranty.

      Take care,
      pavlovsr

  • changsoo lee

    I would like m file 🙂
    I’m learning chemical engineering in Yonsei Univ. Korea.
    Plz help me T T
    ch1ch2ch3@naver.com

  • changsoo lee

    oh, sorry, ch1ch2ch3/naver.com plz T T

  • Elias

    Dear Mr. Lee,
    Is it possible to set the a_matrix as variable, so we ask fsolve to adjust the interaction parameters having a set of experimental data?
    Best regards,
    Elias.

    emonteirofilho\ufsj.edu.br

  • pavlovsr

    @Elias,

    You could!

    %knowing r_matrix, q_matrix, T, x1
    shortUnifac = @(a_matrix) Unifac(r_matrix, q_matrix, a_matrix, T, x1) - gamma_exp_matrix;

    fsolve(shortUnifac, a_matrix);

    Although I wouldn’t use this approach if you don’t have values to compare to… simply because the equations are very non-linear and the value for which you are trying to find is embedded in two exponentials. If you try it let us know!

  • Bobby J

    Hi there.

    I’ve been trying for a few weeks now to get UNIFAC to model a system for me through matlab. So far the best I’ve managed is to get numbers that are in the right ballpark (ie. not out by orders of magnitude from the actual values) but I am having real trouble getting past that.

    It’s for the purpose of a university design project. Please let me know if you can help asap!

    Bobby
    s0828614 / sms.ed.ac.uk

    • Bobby J

      I should also add that what I’m trying to do is model a 7-component system, so if your program is capable of handling that, I would greatly appreciate having a chance to try it out.

      Bobby

    • pavlovsr

      So there could be a couple of issues:

      1. If you’re using fsolve it could be that you need to set your tolerance option when calling it (ie you need to set parameters ‘TolX’ or ‘TolFun’ to > 1E-6… the default).

      2. Make sure you are not using VLE binary interaction parameters for LLE situations and vis versa! They can show similar behavior but have, say an order of magnitude accuracy.

      3. Your system may not be well estimated by UNIFAC. If I remember correctly UNIFAC does not work well for systems with lots of strong bonding forces (ex: Hydrogen bonding) and works best for non-polar systems. Some items you might find interesting are here:
      http://www.aiche.org/uploadedfiles/students/departmentuploads/thermodynamics.pdf
      http://people.clarkson.edu/~wilcox/Design/BYUprop.doc
      http://people.clarkson.edu/~wwilcox/Design/thermodl.htm

      Those can help you decide if UNIFAC is right. In any case I will send my code and add you to the user’s list.

      Good luck,

      Ryan

      • pavlovsr

        A seven component system is possible with this code but it’s not built in… codes like HYSYS and ASPEN I believe just do several binary system estimations for things like this (ie seven components = 49-7 interaction coefficients since there are 7 components and because an interaction coefficient b(i,j) does not necessarily equal b(j,i) ).

      • Bobby J

        Do you mean they take into account that many binary interactions? [ b(i,j)s and b(j,i)s ] or that they only calculate activity coefficients for binary systems and then estimate activity within the mixture? If it’s the latter case that would make modelling this system considerably easier, but I am not sure how I would go about estimating activity in the mixture from activity in binary systems.

        If this is what you meant, could you go into some more detail about how activity coefficients for a system can be built up from those calculated for binary systems?

      • pavlovsr

        So I have never done this myself but something that illustrates how you might go after this is listed in the literature:

        http://onlinelibrary.wiley.com/doi/10.1002/aic.690180208/pdf

        The section you may be most interested in starts at THREE ION SYSTEM part. This may not be what you want to do since the mixing rule that this paper uses to combine the activity coefficients may not be the best description of the system… this is something that you might do at a PhD or Master’s level to get the details right and it would probably be best to consult your thermodynamics professor about this topic. If you do get somewhere I would be interested so let us know!

      • Bobby J

        Thanks for the info, but with deadlines coming up we’ve abandoned trying to get UNIFAC to work through matlab for this system (and found ways around needing activity coefficients – or use UNISIM to simulate non-idealities).

  • Santiago

    I would like the m file for the UNIFAC method for matlab. Thanks in advance!. santiagorodriguez23 / live.com.ar

    • pavlovsr

      Hey Santiago,

      I can give you the code however the code comes with absolutely no warranty and I need to ensure that it is only used for educational purposes. If you can supply me with the name of your learning institution/how you intend to use the code, I can provide you with a copy.

      • Santiago

        Yeah, I study at the University of Los Andes, in Bogotá, Colombia. Right now I’m doing a class (homework) project where we develop the industrial process of Isoamyl acetate manufacturing and scaling, so we first need to find for the quaternary system (acetate+acetid acid+alcohol+water) all the activity coefficients by UNIFAC based in papers and articles that use it. It’ll be great if you could provide me with a copy. as you see is intended to be used for EDUCATIONAL purposes ONLY, cause is just a class project. Thanks for the fast reponse!!

  • Santiago

    Yeah, I study at the University of Los Andes, in Bogotá, Colombia. Right now I’m doing a class (homework) project where we develop the industrial process of Isoamyl acetate manufacturing and scaling, so we first need to find for the quaternary system (acetate+acetid acid+alcohol+water) all the activity coefficients by UNIFAC based in papers and articles that use it. It’ll be great if you could provide me with a copy. as you see is intended to be used for EDUCATIONAL purposes ONLY, cause is just a class project. Thanks for the fast reponse!!

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: