Cron Based Discrete Event Simulator

In my last post I showed how to use the croniter package to find all the cron events that occur between two dates.  I wanted to see if it was suitable for integration into my discrete event simulator, and after a overhaul of the original simulator that's exactly what I've done.  To demonstrate the functionality it's helpful to have a real world scenario to model.

Scenario

A programmer is addicted to energy drinks and every night at 8 pm they put all the empty cans on their desk in the bin.  Some days are more stressful than others so the number of cans taken out is a random number between 0 and 10 included.  The rubbish bin is emptied every Monday morning at 8 am.  The programmer only works from home January through April but the number of cans in the bin needs to be determined between December to May.  The amount of rubbish in the bin needs to be monitored every 30 minutes.

The first step is to define tasks that need to be performed.  They are simply functions that accept a variable describing the simulation environment.

Define Events

The parameters to initialise the simulator include a start time, an end time, and a default location.  The times are made up of a year, month, day of the month, hour, minute, second and an optional location that is a valid entry for the pytz package.  If no location is supplied for the times, the default location is used.  Other parameters specific to the state of the simulation are included as parameters.  Locations are needed as the underlying simulator uses unix timestamps as the timebase.  This prevents ambiguous time due to daylight saving and different time zones.  Essentially a time and date is useless information without a location.

Initialise the simulator

The tasks are scheduled to run via a cron string and a location for the event.  If a location isn't supplied the default simulation location is used.  You can also set a range for the cron events to run.  In the example below the log_data and the empty_rubbish tasks don't have start and end times so they run for the entire simulation.  The take_out_trash task runs for a shorter period of time and has defined start and end times.  The simulation is then run with the results logged to file.

Schedule Events

The results show that on one occasion there are at more than 50 cans in the bin.  This could be used to determine the size of a bin to buy if the current one is too small.  Let's say you have a bin that takes only 40 cans, it wouldn't be too much trouble to calculate the percentage of the time that there are more than 40 cans to be placed in it.  You may decide that it's acceptable for the bin to be overfull for a certain percentage of the time.  It's all about quantifying data from the simulation and using it to make informed decisions instead of just going by your gut.

Simulation Results

This is a work in progress and there are parts I want to improve, and although this scenario is simple I hope I've demonstrated the value of simulating a situation.

Get the Code!

Adding a Sales Profile to a Discrete Event Simulation

If you've been following along, you'll know I've been putting together a discrete event simulator for inventory management.  The goal is to better integrate knowledge of sales history, stock levels, and ordering patterns to make quantitative decisions about what to order and when.  So far I've been using a scenario involving a 24 hour business with a constant sales rate, but that's not realistic.  Therefore the next step is to add a way to generate a customer flow based on a sales profile.

The way I've gone about doing this is to start with a graph of a sales profile.  Let's slightly change the scenario I've been using.  We're now looking at a profile of sales at a fruit shop.  The fruit shop only opens between 8am and 9pm on weekdays, 8am to 5pm on Saturdays, and 9am to 6pm on Sunday. Weekday sales are $3000 a day, and weekend sales are $7500 a day.  This is a total of $30000 per week.  When graphed against seconds passed, the total sales for the week are shown below.

sales profile

In the previous simulations I was modelling the customer flow as a poisson process, where item sales per second was used to generate random times between customer arrivals.  This makes it hard to create a varying sales pattern based on customer flow.  This is easily changed though, you just need to interpret the problem differently.  If for example 150 watermelons are sold per week, that means that on average $200 are spent in the store as a whole between each watermelon purchase.  If we assume that the customer flow of people buying watermelons is similar to those buying fruit in general, which in the absence of better data is a reasonable assumption, a sales pattern can now be generated.

If instead of saying "a customer will buy a watermelon at x seconds into the week" we now say "a customer will buy a watermelon at x dollars into the week", it makes sense to reverse the situation and use the sales as our independent variable.  The way to do that is indicated in the image above.  The evenly spaced horizontal black lines indicate sales increments of $1000.  When projected across to the sales profile and then down to the time axis, a customer flow is generated.  You can see that there won't be any sales when the store is closed as the line can't intersect with these time periods.  Days like Saturday and Sunday that have a high rate of sales intersect a lot of lines and project them into a small time period.

I've written a custom class to perform this interpolation process on a piecewise function that's passed to it.  It allows a dollar value to be converted to a time stamp and vice versa.  The supplied interp function in numpy is almost adequate, but it has some extra constraints placed on it.  Firstly the data supplied needs to have unique timestamps in increasing order, it doesn't make sense to have two different sales values for a single point in time.  The other requirement is that the sales data needs to be increasing, in this case it's fine to have two data points the same, it just means no sales have occurred in this period.

To ensure consistency, I wanted to ensure that when converting from dollars to time, predictable and consistent results were produced.  This may not seem obvious, but if you were to convert the dollar value of 3000 to a timestamp in the above example there are multiple answers.  Any time between closing on the first day and opening on the second day the sales value is $3000.  I decided to make it so the earliest time a particular dollar value is reached would be returned.

It's simply a matter of creating a piecewise function as shown below.  It can be as detailed as you want.  I've only indicated sales at opening and closing times and linearly interpolated between the data points, but you could supply hourly or minute by minute sales data if you wanted.

specification of a sales profile

Using this function, a customer pattern is generated in terms of dollars and then back projected to timestamps.  These values are then used in the simulation.  I haven't changed the ordering requirements from the previous scenario.  You have to order 24 items at a time, the deliveries happen at 1 am, and after the delivery driver has left, you have to have 50 or more items.  So what does the simulation look like now?

stock on hand simulation with sales profile

You can see the plateaus where there are no sales because the store is closed.  The deliveries at 1 am also appear as interruptions to the plateaus.  Also visible are the slow weekdays sales of $3000 over 13 hours.  This can be compared to the weekend sales of $7500 over 9 hours.  This increased rate of sales of approximately 3.6x on a weekend is visible as well.

I'm not happy with my code at the moment so I'll take some time to refactor it and change some of the terminology to generalise some parts.  My main concern at the moment is to create good libraries, the test.py program is sort of a test bed for these.

Get the code!

Monte Carlo Simulation and Discrete Event Simulation

I've been doing some more work on my discrete event simulator, or to be more accurate, the supporting functions that help to analyse the results.  In my previous post I described a basic stock management situation that can be modelled.  I'll continue to use that model here as it's simple to understand and helps illustrate the features I've added.

To bring you up to speed, the simulation models a 24 hour, 7 day a week business that gets deliveries at 1 am and sells items all day.  It's not realistic but it's a starting point.  The main problem with it is the discrete random nature of the events.  This can be seen in the week long simulation below from the last post.  It's possible to identify trends and events, but it offers no insight into how likely these events are to occur.

Original One Week Simulation

Getting data about the likelihood of events occurring can be accomplished by running the simulation for a longer time period and then comparing smaller sections of the simulation.  In my case I ran the simulation for 200 weeks, broke the data up into weeks, and then performed analysis on those multiple trials, just like a Monte Carlo simulation.  These results can be seen below.

The graph illustrates the exact same model, but instead of one trial, the model is essentially run 200 times.  This allows percentiles for the simulation to be generated.  In this case 0, 10, 50, 90, 100.  To put this into more understandable terms, in the graph below, 0% of the simulations predicted stock levels less than the dark blue line, 10% predicted levels less than the green line, 50% predicted levels less than the red line, 90% predicted levels less than the light blue line, and 100% of the simulations predict levels less than the purple line.

A monte carlo analysis of a discrete event simulation

This kind of analysis allows sensible measured conversations about stock levels to take place.  All to often predictions are made that don't allow nuance and can cause problems.  For example, instead of saying "I reckon we'll run out of stock by 1pm"  we can now say "There's a 40% chance we'll run out by 1 pm".  Everyone can take that on board and take appropriate action, it may be a 40% chance is acceptable as there may be some other more urgent matter to take care of.

I'm anthropomorphising the conversation, but in reality I see this as a way for computers to compare and prioritise tasks.

You may have noticed a problem with what I've described.  The original data from the discrete event simulator isn't a regularly sampled data set.  A customer may come in 10 minutes after the last one and then you won't see another for an hour.  To perform the Monte Carlo simulation the data needs to be resampled with a regular timebase.  In my case I've choose to sample the data every 5 minutes.

This is done with a kind of zero order hold.  The process to determine the value at a specific time is straight forward, if an event occurs at this time the new data point will be the value after the event.  If there is no event, the data point will be the value after the last event.  If there is no previous event the first one is used.

The image below illustrates the initial data set of (1, 1), (5, 2), (10, 3), (11, 4), (12, 5), (20 , 6).  These are marked with blue crosses.  The red line shows how this data is resampled for different points.  The black bars show the result of resampling the data at -10, 1, 2, 10 and 30.  The results are 1, 1, 1, 3, and 6 respectively.  I've chosen to resample the data at 5 random points, but if performed at regular intervals the information is ready to be fed to the Monte Carlo analyser.

Resampling a data set

To give a usage example I've folded the entire simulation into a period of one day.  That means there are 1400 days sampled every 5 minutes simulated and compared in the image below.  If I had concerns about the stock levels at 8:50 pm (sample 250), I could look at this and say that I'm very confident that there will be stock, if fact there should be about 30 to 50 items available.

The simulation folded into a single day

In the simulation below I've increased the sales rate so a different question can be asked.  If customers complain about out of stocks you can look at this and say that supplies will run out somewhere between 7:30 am and 3:50 pm (90 and 190), but most likely at 10:50 am (130).  This then allows you to take action, such as increasing the order size.

Increase rate of purchase indicates likely stock exhaustion points

It's all a bit simplistic at the moment, but it's getting there.  What I'm really aiming to show is the complex interaction of seemingly unrelated variables.

Get the code

.

A Basic Discrete Event Simulation

Let me paint a picture.  You run a service station (or gas station for our American friends), and you sell bread.  It's a 24 hour business with a uniform flow of customers throughout the week and you want to investigate the effectiveness of the ordering routines you use.  How do you do it?  The answer is Discrete Event Simulation.

I'm ashamed to say I've never heard of it, and once again came up with the idea independently  (Another case of "Wouldn't it be cool if ........  *Googles*  Oh, this already exists, is well developed, and people make a living from it").  I'm familiar with the simulations used in designing electronics, but had never experienced this type.  Basically a queue of events with timestamps is created, the next event is removed and processed and this is repeated until a certain time or all the events have been processed.

The reason I'm looking into this is because I have something I want to simulate and had originally planned to use the Simpy simulation package for python, but as I started working it became apparent that it was overkill for what I needed.  It's great for simulating resources that interact with each other, and although it looked like my problem needed this it really didn't.  The ordering process and customer process don't interact, even though it looks like they do.  I'm referring to the time that events occur.  Customers arrive indepenently of each other and of when deliveries happen.  They are linked by the stock levels in the store but this is handled by a stock on hand variable (My basic scenario above could probably be handled by Simpy, but I plan for it to get complicated).

After much procrastination I decided to write my own, but as it turns out someone has done the hard work for me.

User GaretJax  from http://codereview.stackexchange.com/questions/3670/discrete-event-simulation-with-variable-intervals already had a basic framework for this problem.  So I ran with it.  (Check github for  DiscreteEventSimulator.py)

The first thing to do is describe the processes at play.  The customer for example is described below.  They walk in and if bread is on the shelves, they buy a loaf and the stock on hand is decremented.  The time until the next customer is described by an exponential distribution to simulate a poisson process.  The final thing to do is schedule the next customer arrival.  The sales rate is 0.0004 loaves per second, or more intuitively, one loaf every 2500 seconds or around 40 minutes.

Customer Process

The ordering process is very simple.  The bread delivery truck turns up and the staff member counts the stock the store currently has, and then works out how many loaves are needed to bring the stock levels back to 50 loaves.  But there is a small problem.  Loaves of bread only come in crates of 24.  So they need to round up to the nearest multiple of 24 and take that from the driver.  This is then added to the stock on hand.

Order Process

All that's needed now is to start the simulation.  7 events are created for the bread deliveries, one for each day of the week at 1 am.  An initial customer is also scheduled.  Only one is needed as each customer generates the next customer event.  I've then set the simulation to run for a week.

Simulation configuration

The results below starting from Monday are interesting.  We can see the random nature of the purchases, but the overall trends are visible.  We can see that on Monday, Tuesday, Wednesday, Thursday, Friday, and Sunday, only one crate a day was needed, but on Saturday two crates were needed.  It's also possible that they're holding too much stock as for the first 3 days of the week the stock barely falls below 40 loaves.  (Customers have been complaining about stale bread, maybe this is the reason.)

Simulation Results

This is only a quick overview of the process with some very simplistic assumptions and modelling, but I've got a framework to use now and I have thoughts on how to make this more realistic (and complicated).  This problem has been kicking my arse all week,  I only managed to pull it all together in the last couple of hours.  Very happy at the moment

Get the code!

LED Driver Boost Converter

So, it's time for an update on the LED driver I'm designing.  I think I have the final design of the circuit and have most of the components selected.  I haven't bothered selecting actual parts for non critical components yet.  The design is mostly locked in, a value may change here or there but it's pretty much final.  The input filter added to remove conducted EMI may need some improvement as there may be the possibility of ringing.

I'll admit that this has taken longer than I would have liked, but I took it as an opportunity to learn how all the components of a boost converter interact and change the overall operation.  I didn't want to blindly follow a reference design and not know what I was doing.  I still have a lot to learn, particularly about compensation and stability, but it's been fun.

Boost Converter LED Driver - with PWM and Analog Dimming

The plot below shows the circuit starting up.  The PWM dimming function is in use with the frequency being swept from 0 to 3kHz.  The output slew rate limiting circuit does its job and brings the output voltage up gradually without causing too much current to flow through the main inductor.

LED Driver Sweeping the PWM Frequency from 0 to 3 kHz over 60 ms

In the plot below you can see a cycle of the converter operating.  The LED currentin dark blue slightly overshoots its target but comes into regulation rather quickly.  The main inductor current is shown in light purple and can be seen to have ripple, but nothing that wasn't anticipated.  The main supply current is shown in dark purple and has no noticeable high frequency component, but there is some ringing,  I don't want that.  So this area needs some attention.  If anyone has experience with this I'd be glad to take suggestions.

I've simulated the snot out of this thing, but that can only get you so far.  So it's time to actually do a board and see what happens.

Main Inductor Current, Supply Current and LED Current

You can grab the associated files from the GitHub Repository

Boost Converter Output Capacitor Selection

I was a day late with this post as I was trying to select capacitors for my boost converter.  I knew it wouldn't be easy, but I didn't think It would be this hard.  The loose specification I'm trying to achieve is listed in this post, but basically I'm trying to convert 12 to 60 Volts driving a 300 mA load.

I started by looking at different converter topologies, and decided to go with a boost converter due to what parts I could get (I didn't want to hand wind transformers - walk before I can run).  I'm reasonably sure it's not the optimum topology, but it should be doable with off the shelf parts.  I'm at the point where I've got a stable design that does what I want, and am now trawling through suppliers to find actual parts.  Things were going great,  got an inductor, got a switching Mosfet, got a Schottky diode, hit a wall when it came to the output capacitor.

The problem is that the capacitor needs to be rated for 60 Volts, and to minimise ripple needs to be low ESR, and has to be around 5 uF.  That's like trying to find pixie dust.  Let's look at the equations that determine what capacitor I'm looking for and see if we can't tweak some things, keeping in mind that the purpose of the output capacitor is to minimise ripple in the output voltage.

Levels of Output Voltage Ripple in a Boost Converter

The output voltage ripple consists of two parts, the first part occurs when the capacitor discharges while suppling the load while the Mosfet is on and charging the inductor.  This means it's related to the duty cycle, frequency of operation, load current and the capacitor size.  I can't change the load current or the duty cycle, but I did decide to up the frequency of operation from 400 kHz to 600 kHz.  This means I can get away with a smaller capacitor.

The second part of the voltage ripple is caused by the ESR of the capacitor. When current flows into and out of the capacitor through the ESR during its charge discharge cycle, a ripple is created.  Reducing this ripple can be done in two ways, decrease the ripple current in the inductor, or get a low ESR capacitor.  I've increased the inductor size to reduce it's current ripple, but the best option is to select a low ESR capacitor.  Ceramic capacitors have such a low ESR they are essentially considered to have none, and for this reason have become a favourite in DC-DC converter designs.  When going through all my searches and considering specs against cost, they're pretyy much the best option.  They have one drawback though, when biased with a DC voltage their capacitance is reduced.  In some cases, significantly.  I'll let EEVblog explain.

So finding the capacitance I need isn't easy, particularly when I need to find parts rated to 60 Volts.  In the end I decided to go with 4 x 12101C225KAZ2A from AVX.  They're X7R 2.2 uF capacitors rated to 100 V each and when biased to 60 V I should get at least 550nF out of them, giving a total of 2.2 uF.

I need to investigate how the dynamic resistance of the load interacts with all of this, but in the end I think everything should work.  I've upped the frequency of operation, and slightly compromised on my acceptable level of ripple.  After all, isn't engineering is the art of intelligent compromise?

Boost Converter Design

.

Boost Converter Startup

My ongoing series on designing a boost converter continues.  Today's post will be short one.  What happens in the converter during startup?  The period between when you first apply power and when the controller starts to do it's thing.

Using the schematic further down the page as a reference, the startup waveforms below can be described.  It's also helpful to know the period of an LC oscillator as well.

steady state period of an LC oscillator

1. During startup the power mosfet (M2) is off, and can be ignored in this analysis.  The current in inductor L1 (36 uH) is zero and the voltage across the output capacitor bank C5, C7, C8, C9 (9.9 uF) is equal to zero.  When the 12 volt supply is first applied, the inductor L1 and the capacitor bank start to resonate.  This causes the initial peak.  The initial inductor current (purple) should resonate around zero amps.  Knowing that the peak voltage deviation will be 12 volts (i.e. the voltage across the capacitor will have a range of 12 ± 12 volts) it can be shown that the current through the inductor L1 will be equal to 12*sqrt(C/L).  In this case that equals 6.29 Amps.  We get close with around 5.8 Amps.  The forward voltage of the diode means that the supply is more like 11.4 Volts and there are also some series resistances that will cause the current to be lower.  That equation is only valid for steady state analysis as well.
2. In theory the voltage across the capacitor bank (red) should peak at double the input voltage, as an undamped LC circuit will swing to double the magnitude of a step input to give 24 Volts.  In this case however, the forward voltage of the diode and small resistance of the inductor has reduced this to the 22.4 ish volts observed in the simulation.  At this point we have approximately one half period of oscillation.  Although meant for a steady state conditions, this can be checked against the formula above.  t = 2*pi*sqrt(LC) = 118us.  Half a period is approximately 60 us, pretty close.  Another effect starts to happen at this point as well, for the first time the inductor current goes to zero and tries to reverse.  The diode doesn't allow this to happen, leaving the output capacitors charged at around 22.4 volts.  The inductor starts to immediately resonate with the parasitic capacitance of the mosfet.
3. The damped oscillations have a period of approximately 1.05 us.  Using the equation above it can be shown that the parasitic capacitance is 775 pF.
4. The current ripple can also be seen in the inductor.

(Click for larger images)

Startup waveforms of a boost converter

Boost converter

Get the associated files here.

How Can Current Flow Backwards Through The Inductor of a Boost Converter?

I've been keeping myself occupied by (slowly) designing a boost converter for an LED light panel, and during one of my simulations I noticed something weird, there was current flowing backwards through the inductor.  I couldn't seem to find any reference to this anywhere so I thought I'd write up what I think is happening.

First I'll give a quick explanation of how a boost converter is meant to work using the schematic below.  When the power mosfet M2 conducts, the input voltage is placed across the inductor L causing the current through it to rise.  When M2 is switched off, the current through the inductor now flows through D21 to fill the capacitor C5 and supply the load.  I won't go into the maths, but it's trivial to calculate the voltage gain of the converter.  It's important to note that current can't flow backwards through the inductor at any point.  For this to happen, current has to either go backwards though the diode D21 (which happens a tiny bit, but isn't what we're seeing) or it has to come out of the mosfet M2, which doesn't make much sense either.

Boost converter schematic

The simulation that initially caused me to look into this is below.  It occurs during startup, when the output capacitor isn't fully charged and when the controller is in soft start mode.  At this point the circuit is operating in discontinuous mode, meaning that the inductor current will fall to zero during each switching cycle.  Nearly every textbook will illustrate this with a triangular waveform with a flat section in between triangles to indicate the zero current section.  This isn't what happens though.  I'll step through the stages of operation in the image below.

1. The mosfet is turned on (gate voltage is the green waveform) and the current in the inductor starts to ramp up (light blue).  The voltage at the point between the inductor and mosfet drops to zero (red) and there isn't any current flowing through the diode (dark blue).
2. The mosfet is turned off.  The current flowing through the inductor starts to ramp down.  The voltage at the point between the inductor and mosfet increases to near the output voltage. It can be seen that the inductor current goes through the diode as the dark and light blue plots are the same in this section.
3. This is the point where traditional explanations fail.  The mosfet is still off and the diode current has reached zero.  The voltage between the mosfet and inductor starts to fall as you would expect.
4. The inductor current keeps ramping down though, and in this case reaches around 200 mA in the reverse direction.  The only place this can come from is the mosfet.

Boost Converter Simulation

So what's going on here?  I contend that this can be explained by taking the parasitic capacitance of the mosfet and diode into account.  I'll leave exactly what capacitance causes this up to someone more familiar with semiconductor physics.  To demonstrate this and remove other variables I've put together a simplified model below.

The input voltage source remains. As the output capacitor is made increasingly large it starts to act like a voltage source, therefore it can be replaced with a voltage source set to the output voltage we are interested in, and a steady state analysis can be performed.  The mosfet has been replaced by a voltage controlled switch with a parallel capacitance, in this case 500 pF, to represent the parasitic capacitance.  A reverse biased diode has also been placed in parallel to represent the mosfet body diode.  The Schottky diode connected to the output has been replaced with an ideal diode with an appropriate forward voltage.

Simplified Model

Note: I apologise for the graphs, you may need to click on them to see things more clearly.  The voltage and current scales don't have a common zero either.

Running this simulation gives an almost identical plot to the one above, but we can now see what is happening more clearly.  The inductor ramps up as before (dark blue), it then ramps down to zero amps and the output diode stops conducting (red), but we can now see that the parasitic capacitance is still charged (green), it starts to discharge, but can only do this via the input inductor.  This is the reverse current through the inductor.  What's happening is that the parasitic capacitance is resonating with the main inductor around the zero volt point.  In the simulation below the resonance occurs at a frequency that makes it look like it's part of normal operation.  The period of resonance for an LC circuit is 2*pi*sqrt(LC), in this case that's equal to approximately 0.5 us.

It's hard to see but you may also notice a tiny delay between the time the mosfet is turned off and when the diode starts to conduct.  This is the first bit of current from the inductor charging the parasitic capacitance.  Once this capacitance is charged and the voltage reaches a high enough point, the diode will then start to conduct.

Simplified simulation

All of this can be made clearer by shortening the on period as seen below.  Ideally the time spent with a negative inductor current should be around half a period of the resonance 0.25 us, which is close, but there is another effect at play that can be see in the green trace.  The voltage between the inductor and mosfet hits a limit of -0.7 volts.  At this point the body diode starts to conduct, effectively dampening the resonance.

Simplified simulation - reduced on time

If we were to reduce the parasitic capacitance by a factor of 5, the resonance period would be reduced by a factor of sqrt(5), giving 0.22 us.  Which is once again close to the simulation below.  This time the resonance is clearly seen in the voltage between the inductor and mosfet.  It can also be seen where the bottom of a typical resonance waveform has been clamped by the body diode.

Simplified simulation - lower capacitance

This is another one of those theoretical cases.  I don't know if this will really happen when a real device is built.  I don't even really know if it's a bad thing, I just had to put an explanation up in case someone else has the same issue.

You can get the related files here.

Edit:  On reflection, it is a bad thing.  High frequency noise is a bad thing in any circuit and leads to EMI problems.  I've since found out that googling "Boost Converter Switch Ringing" will yield more information about this issue.